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Non-perturbative Methods in Modal Field Theory 

Abstract 

Several issues in the modal approach to quantum field theory are discussed. Within 
the formalism of spherical field theory, differential renormalization is presented and 
shown to result in a finite number of renormalization parameters. Computations of 
the massless Thirring model in 1+1 dimensions are presented using this approach. 

Diagonalization techniques in periodic field theory are demonstrated. Issues of 
very large Hilbert spaces are considered and several approaches are presented. The 
quasi sparse eigenvector (QSE) approach takes advantage of the relatively small num- 
ber of basis states that typically contribute significantly to any particular eigenvector. 
Stochastic correction methods use Monte Carlo calculations to calculate higher order 
corrections to the quasi sparse result. 

The quasi sparse eigenvector method and stochastic error correction are applied 
to the Hubbard model. With j — 4, the shift in the ground energy below the U = 
value is found to within 1% for the 8x8 Hubbard model with || filling. 
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Chapter 1 
Introduction 



I first encountered quantum field theory as an undergraduate. I was drawn by the 
beauty of its symmetric and seemingly simple equations. But I was also enticed by 
the opacity of those same equations which so stubbornly resist calculations. 

Free theories can be calculated exactly and others can be calculated in perturba- 
tion theory, but one quickly runs into problems of infinities. Even after regularization 
and renormalization the perturbation series is at best asymptotic and for many phys- 
ical systems is virtually useless. 

Lattice regularization is another approach which avoids the constraints on coupling 
constant. From a constructivist viewpoint, it is reassuring that the field theory can be 
put on a lattice and a finite answer extracted. The fermion derivative term presents 
problems in discretized space. Although these can be handled, the complications thus 
engendered invite a new approach. 

Thus, when I was introduced to spherical field theory, I was initially attracted by 
two main features. The first was that although the answer would still be expressed as 
an infinite sum, just as in perturbation theory, the series would in principle converge. 
The second was that, since space is still treated as a continuous variable, fermions 
could be treated in a naive manner. 

On the other hand, since, space was still continuous, spherical field theory faced 
the problem of renormalization. Because the Hamiltonian and the natural regulators 
are functions of "t", the radial coordinate, and thus, it is conceivable that arbitrary 
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functions might be required to renormalize a theory. Diagrammatic renormalization 
is only useful for super-renormalizable theories with a finite set of such diagrams. The 
second chapter, "Renormalization in spherical field theory", addresses this problem. 
It shows that a small set of local counterterms is sufficient, in general, to remove all 
ultraviolet divergences in a manner such that the renormalized theory is finite and 
translationally invariant. 

The Thirring model in 1+1 dimensions has a four fermion interaction and is not 
super-renormalizable. The next chapter, "The massless Thirring model in spherical 
field theory" , serves as concrete test of the regularization scheme. It also served as 
a laboratory to test the efficacy of different techniques for handling fermions. The 
Hilbert space of the system was small enough to fit in memory and a direct Runge- 
Kutta approach was used to integrate the equations of motion. 

A Monte Carlo integration was also attempted at this time but the results were 
mediocre and were not included. Some difficulties of the spherical approach became 
apparent during this investigation. The time dependence of the Hamiltonian made 
small steps necessary near t — but at large t, the Hamiltonian is small and a long 
time is necessary for the state to evolve to the ground state. This problem is mainly 
technical and was solved with adaptive time steps. 

The time dependence in the regularator also created a moving target for the 
number of modes necessary for the calculation. At small t, where Mt was small 
only a few modes would be necessary but at larger t, more were required. The 
majority of computer time was spent calculating the time evolution for small t. The 
requirement to include the extra modes for their effect at large t resulted in a waste of 
computer resources. Finally, the time dependence of the Hamiltonian prevented the 
precalculation of certain constants and other optimizations that a time independent 
H would have allowed. 

In quantum mechanics, the greatest optimization for calculating time evolution is 
expression of the system in terms of energy eigenstates. The next chapter tackles this 
problem head-on for 4 theory in 1+1 dimensions. Space is now a periodic box of 
length 2L. The advantage of a time independent H has been gained at the cost of a 
new parameter, L which must also be taken to infinity before we can use our results. 
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The work in this chapter was exciting to me for two reasons. One was the potential 
to work directly in Minkowski space. Even more exciting was the potential for explicit 
examination of particular eigenstates of the system. Figure 8 in chapter 5 is a perfect 
example of the types of work that could be done within this approach. The eigenstates 
are tracked across a phase transition and the symmetry relationships are exposed. 

At this point we were ready to try our new approaches on systems in 2+1 or 3+1 
dimensions. The exponentially higher dimensionality of the Hilbert space became the 
dominant issue we faced. I considered simply waiting. It is a "rule" in the computer 
industry that computer power doubles about every 2 years. If the rate sped up a 
little, in thirty years we might be able to attack some 2+1 dimensional systems and 
in another thirty years we could try problems in 3+1 dimensions. 

The following question then presents itself. "Are these systems in 2+1 dimensions 
so complicated that they cannot be described without reference to 10 18 states or 
more, or after the diagonalization was completed would we find that the results could 
be described in a simpler way?" In particular, it is possible that any particular 
eigenstate of the Hilbert space may require only a relatively small number of basis 
states to accurately describe it. 

Chapter 5 argues for the prevalence of this condition, known as "quasi-sparsity" . 
A careful analysis shows that this is as much a statement about the types of bases 
we are likely to use as it is a statement about the Hamiltonians we will encounter. In 
our work we use either a Fock state basis or a position state basis. Different problems 
turn out to fit into the quasi-sparse model to different extents. As pointed out in 
chapter 8, the Hubbard model with a Fock state basis seems particularly ill-suited 
for this approach. 

The power of lattice field theory comes from its ability to tap into Monte Carlo 
computational methods. The dimensionalities of the relevant spaces are even larger 
than those considered in this thesis but because of the smoothness of contributions 
as a function on the configuration space, it is possible to do a good job of sampling 
the important configurations. I am convinced that efficient use of Monte Carlo is 
essential to solving most physically relevant quantum field theory problems. 

In lattice field theory, Monte Carlo is restricted to imaginary time calculations 
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and cannot handle unquenched fermions well. A successful integration of Monte 
Carlo computation into the approach described in chapter 4 could address these 
limitations and give visibility to more than one state in a symmetry sector. Two 
potential methods of doing this integration are presented in chapter 6. 

The thesis ends with an attempt to apply the results of chapters 5 and 6 to 
the Hubbard model. It was expected that the complexity of its ground state would 
present a challenge for our methods and we were not disappointed. Adjustments to the 
stochastic methods of chapter 6 are presented in chapter 8 and a reasonable estimate 
of the ground state energy of the 8x8 Hubbard model is extracted. Other questions 
such as the value of particular correlators in the ground state could also be computed 
with this framework. Other information, such as the makeup of excited states will 
require further modifications. Some potential directions for further improvement are 
mentioned in the chapter. 

One area where our current state of the art seems particularly lacking is in the 
sampling of configurations for Monte Carlo. We run trajectories so the choice at 
each step has no global knowledge of its path. While I was able to find a good 
"distance" based guiding scheme for the Hubbard model it does not take the energy 
of the states into account. It may be that the Hilbert space of a fermion system 
does not support the same notion of close paths that the lattice field theory boson 
computations use. But if it does, a sampling scheme which uses it will almost certainly 
be an improvement. 

The order of chapters in this thesis corresponds with the chronological order of 
the papers they contain. They also mirror the logical progression of my approach to 
quantum field theory. In trying to do computations on more difficult systems I have 
been driven to adopt features of the lattice field theory approach. The goal for the 
future will be to keep some of the advantages of modal field theory while learning 
from the sampling techniques of lattice Monte Carlo. 



Chapter 2 

Renormalization in spherical field 
theory i 



2.1 Introduction 

Spherical field theory is a non-perturbative method which uses the spherical partial 
wave expansion to reduce a general <i-dimensional Euclidean field theory into a set of 
coupled radial systems (|3|| 36, [|]). High spin partial waves correspond with large 



tangential momenta and can be neglected if the theory is properly renormalized. The 
remaining system can then be converted into differential equations and solved using 
standard numerical methods. 4 theory in two dimensions was considered in In 



that case there was only one divergent diagram, and it could be completely removed 
by normal ordering. In general any super-renormalizable theory can be renormalized 
by removing the divergent parts of divergent diagrams. Using a high-spin cutoff J max 
and discarding partial waves with spin greater than J max , we simply compute the 
relevant counterterms using spherical Feynman rules. 

The J max cutoff scheme however is not translationally invariant. It preserves ro- 
tational invariance but regulates ultraviolet processes differently depending on radial 
distance. In the two-dimensional 4 example it was found that the mass counterterm 



X D. Lee, N. Salwen, Phys. Lett. B460 (1999) 107. 



5 



Chapter 2: Renormalization in spherical field theory 



6 



had the form 

oc 2 (t) [K (rf)I (fJ,t) + 2E b= i,j3M4W] , (2-1) 

where J„, are n th -order modified Bessel functions of the first and second kinds, \i 
is the bare mass, and t is the magnitude of t. As J max — > oo, we find 

C c ,. ex 2 (t) [log(^) + 0(J m l x )] . (2.2) 

Our regularization scheme varies with t, and we see that the counterterm also depends 
on t. The physically relevant issue, however, is whether or not the renormalized theory 
is independent of t. In this case the answer is yes. Any t dependence in renormalized 
amplitudes is suppressed by powers of J^l x , and translational invariance becomes 
— > oo. 

We now consider general renormalizable theories, in particular those which are not 
super-renormalizable. In this case the number of divergent diagrams is infinite. Since 
we are primarily interested in non-perturbative phenomena, a diagram by diagram 
subtraction method is not useful. In the same manner strictly perturbative methods 
such as dimensional regularization are not relevant either. Our interest is in non- 
perturbative renormalization, where coefficients of renormalization counterterms are 
determined by non-perturbative computations.0 In this paper we analyze the general 
theory of non-perturbative renormalization in the spherical field formalism. In the 
course of our analysis we answer the following three questions: (i) Can ultraviolet 
divergences be cancelled by a finite number of local counterterms? (ii) Can the 
renormalized theory be made translationally invariant? (iii) What is the general 
form of the counterterms? 

The organization of the paper is as follows. We begin with a discussion of differ- 
ential renormalization, a regularization-independent method which will allow us to 
construct local counterterms. Next we describe a regularization procedure which is 
convenient for spherical field theory. In the large radius limit t —>■ oo our regular- 
ization procedure (which we call angle smearing) is anisotropic but locally invariant 

2 We should mention that Pauli-Villars regularization is compatible with non- 
perturbative renormalization. However this introduces additional unphysical degrees 
of freedom and tends to be computationally inefficient. 
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under translations. For general t we expand in powers of t~ l to generate the general 
form of the counterterms. We conclude with two examples of one-loop divergent di- 
agrams. We show by direct calculation that the predicted counterterms render these 
processes finite and translationally invariant. 



2.2 Differential renormalization 

Differential renormalization is the coordinate space version of the BPHZ method.^ 
It is framed entirely in coordinate space, and renormalized amplitudes can be defined 
as distributions without reference to any specific regularization procedure. Differential 
renormalization was introduced in [fT3|| , and a systematic analysis of differential renor- 
malization to all orders in perturbation theory using Bogoliubov's recursion formula 
was first described in The usual implementation of differential renormalization 



is carried out using singular Poisson equations and their explicit solutions. In our 
discussion, however, we find it more convenient to operate directly on the distribu- 
tions.0 We describe the details of our approach in the following. We should stress 
that the two approaches are equivalent, differing only at the level of formalism. 

We assume that we are working with a renormalizable theory. For indices • • • ij 
let us define 

t n,-ij = t i n i2 ---t i ', (2.3) 
V h ,.. ij =V h V i2 ---V ij . (2.4) 

Let fit) be a smooth test function, and let I(t — F; fi 2 ) be a smooth function with 
support on a region of scale ^T 1 . We define Si [f] (t) as I(t — r ; /i 2 ) multiplied by the 
jth orc } er term in the Taylor series of f(t) about the point v. Inserting delta functions, 
3 Paraphrase of private communication with Jose Latorre. 

4 Our approach is similar to the natural renormalization scheme described in [[54] 



In contrast with [|54]], however, we do not a priori specify the finite parts of amplitudes. 
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we have 

(2.5) 

= I(f-?;A* 2 ) £ ^^JV':'V;: ... ;/ r (/' • :•)/(:). 

For the purposes of this discussion we will require 

lit- F; li 2 ) = 1+ N {t- ?) as F -> £ (2.6) 

where iV is some positive integer greater than the superficial degree of divergence of 
any subdiagramP] in the theory we are considering. For any renormahzable theory 
N > 2 will suffice. In our formalism, I(t — v\ /i 2 ) determines how finite parts of 
renormalized amplitudes are assigned, and /i is the renormalization mass scale. 

We now consider a particular diagram, G, with n vertices. We define K(ti, ■ • • t n ) 
to be the kernel of the amputated diagram, i.e., the value of the diagram with vertices 
fixed at points t\, ■ ■ - t n . The amplitude is obtained by integrating K(t\, ■ ■ - t n ) with 
respect to all internal vertices. We will regard K as a distribution acting on n smooth 
test functions fx,---f n - (For external vertices containing more than one external line 
and/or derivatives, f e xt(t e xt) should be regarded as a product of test functions, with 
possible derivatives, at t ext .) 

K:f h ---f n ^J d% ■ ■ ■ d% K(t h ■ ■ ■ Qhih) ■ ■ ■ f n (t n ). (2.7) 

Let us assume that our diagram is primitively divergent with superficial degree of 
divergence j. We now define another distribution TqK, which extracts the divergent 
part of K. We start with the case when G has more than one vertex. Let us define 
T G K:f ir .-f n ^ 

E Jd% ■ ■ ■ d% K(t h ■ ■ ■ t n )Sp Mh) ■■■Sf f n (t n ), (2.8) 

JlH yjn<3 

5 In our discussion a subdiagram is a subset of vertices together with all lines 
contained in those vertices. 
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where t ave = -(t% + ■ ■ ■ + t n ). We note that the subtracted distribution K — T G K is 
finite and well-defined for all fx, ■ ■ ■ f n . Let us define 



(2.9) 



Jd% ■ ■ ■ d% ^ ( *i±_±4 _ t)K{t x , ■ ■ ■ t r . 



We can then rewrite TqK : fi, - ■ ■ f n 



n 

fc=l,---n 



Jk ■ 



E E 

JlH \-jn<3 ix^\,i2,V 

*1 ,n ' "tjn ,n 



jd'ir;. : i?: (t) / d% ■ ■ ■ d 4 z n 

U k=1 ,..n VL,.... J\t- Z k )) h{z X ) ■ ■ ■ f n (z n 



(2.10) 



The delta functions make this kernel completely local. We can read off the corre- 
sponding counterterm interaction by functional differentiation with respect to each of 
the component functions of f ex tit ext ) for the external vertices and setting fi n t{t int ) = 1 
for the internal vertices. We now turn to the case when G has only one vertex. For 
this case we set TqK = K, which is equivalent to normal ordering the interactions 
in our theory. In this case K is itself local and therefore TqK and our counterterm 
interaction are again local. 

We now extend the definition of Tq in (|2.10|) to include the case of subdiagrams. 
Let G be a general 1PI diagram, and let G' be a renormalization partQ of G with 
superficial degree of divergence f. For notational convenience we will label the vertices 
of G so that the first n' vertices lie in G'. If G' has only one vertex then again we 
normal order the interaction. Otherwise we define TqiK : fi, ■ ■ ■ f n — > 



■ d A t n K(ti, ■ ■ ■ t r 



Sf 1 fl(h) ■ ■ ■ Sl n f n '{t n i 

L ave L ave 

'fn'+l(tn'+l) " ' ' fn(tn) 



(2.11) 



where t ave = h{ti + • • • + t n /)Q This definition can be used recursively to define 

products of TqiTqi^ for disjoint subdiagrams G[ fi G' 2 = or nested subdiagrams 

G[ D G' 2 . For the case of nested subdiagrams we always order the product so that 

larger diagrams are on the left. 

6 A renormalization part is a 1PI subdiagram with degree of divergence > 0. 
7 After applying T G /, we regard G' as being contracted to single vertex at t ave . 
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It is straightforward to show that the T operation acts as the identity on local 
interactions and thus treats overlapping divergences in the same manner as BPHZ. 



Following the standard BPHZ procedure ([§], |64|, gtj), we can write Bogoliubov's R 
operation using Zimmerman's forest formula, 

^EnK), (2.12) 

F 7£F 

where F ranges over all forestsQ of G, and 7 ranges over all renormalization parts of 
F. In the product we have again ordered nested subdiagrams so that larger diagrams 
are on the left. Let j be the superficial degree of divergence of G. The renormalized 
kernel, RK, is given by 

RK = RK for j < 

2.13 

RK = (1 - T G )RK for j > 0. 

Our final result is that all required counterterms are local, and the form of the coun- 
terterms is 

£c, = £^v 4 « F A(£)A(<j>(?), V^(t)), (2.14) 

where the sum is over operators of renormalizable type. For the case of gauge theories, 
our renormalization procedure is supplemented by the additional requirement that 
the renormalized amplitudes satisfy Ward identities.^ If our regularization procedure 
breaks gauge invariance these identities are not automatic and the required local 
counterterms will in general be any operators of renormalizable type (not merely 
gauge-invariant operators). This is, however, a separate discussion, and the details 
of implementing Ward identity constraints will be discussed in future work. 



2.3 Regularization by angle smearing 

In this section we determine the functional form of the coefficients Fj^it) in (|2.14j ). 

To make the discussion concrete, we will illustrate using the example of massless 4 

8 A forest is any set of non-overlapping renormalization parts. 
9 See [0, |TT| for a discussion of gauge theories using the method of differential 
renormalization. 
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theory in four dimensions 

£ = l0V 2 0-A0 4 + £ c . i .. (2.15) 

From Q2.14Q C c ,t. is given by 

i^(t)0 2 (t) + Ey^v*(*5Vi^(*>^(*) + *>(*V(*)- (2-16) 

Let G(t,t) be the free two-point correlator. We will use a regularization scheme 
which preserves rotational invariance and is convenient for spherical field theory, but 
one which breaks translational invariance. We regulate the short distance behavior 
of G by smearing the endpoints over a radius t spherical shell within a conical region 
R M 2(t), where R M 2(i) is the set of vectors u such that the angle between t and u is 
between —-hri and Xr (see Figure |2"TT| ). The result is a regulated correlator 




Figure 2.1: Sketch of the angle-smearing region (three-dimensional rendering) 

G M 4tJ) = 1 arr L n ^ d 3 ud 3 u'G(tu,t'u'). (2.17) 

u'eR M 2(?) 

We recall that our renormalized theory is determined by the translationally invariant 
function I(t — t ave ; /i 2 ) described in the previous section. Even though our regular- 
ization scheme breaks translational invariance, the renormalized theory nevertheless 
remains invariant. 
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As the radius t increases the curvature of the angle-smearing region becomes 
negligible. In the limit t — > oo the region becomes a flat three-dimensional ball with 
radius lying in the plane perpendicular to the radial vector. In this limit our 
regularization is invariant under local transformations and the counterterms converge 
to constants independent of t, 

limF^il = M 2 $(iJ) (2.19) 
Imi^ffl (2.20) 

We have chosen our coefficients to be dimensionless. Although our regularization 
scheme is invariant under rotations about the origin, the radial vector has a special 
orientation which is normal to our three-dimensional ball. Our regularization scheme 
is therefore not isotropic. The result (as should be familiar from studies of anisotropic 
lattices) is that the coefficient of the kinetic term has two independent components 

Starting with the t — > oo result at lowest order, we now expand our coefficient 
functions in powers of 

pij (f) _ J3,(P) (J?_\ _j_ J_ iJ.(l) , _J_ iJ.(2) /J^\ _i (2 22) 

= M^(i) + f c£>(&) + + • • • (2.23) 

= + + + • • • • ( 2 - 24 ) 

For the moment let us assume t > A^ 1 for 

A = mlM x ~\ (2.25) 

for some fixed mass m and constant z such that < z < \. In this region our 
dimensionless expansion parameter is bounded by ^ and therefore diminishes 
uniformly as M — > oo. 

In general the -j^j dependence in the functions will contain analytic terms as 
/i 2 — > as well as logarithmically divergent terms. There are, however, no inverse 
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powers of j^. These would indicate severe infrared divergences not present in the 
processes we are considering, as can be deduced from the long distance behavior of 
the integral in (CT-H With this we can neglect terms which vanish as M — ► oo, 

*>(t) = M^(fy + (2.26) 
Fv<t>v<t>$) = c V4>v<i>{'M^) (2.27) 
= (2-28) 

Since our regularization scheme is invariant under M — > —M, we have also omitted 
the term proportional to c}} which is odd in M. 

We now consider what occurs in the small region near the origin, t < A -1 . For 
the theory we are considering (and in fact for any renormalizable theory) the highest 
ultraviolet divergence possible is quadratic.^ In the limit M — > oo we deduce that 
each Fa scales no greater than 0(M 2 ). On the other hand the volume of the region 
t < A -1 diminishes as 0(M 4z ~ 4 ). Thus the total contribution from the region t < A -1 
scales as 0(M 4z ~ 2 ) and can be entirely neglected. 

To summarize our results, the counterterm Lagrange density has the form 

4VV0(t)) 2 + 4%$ ■ V<^1) 2 + (^cj? + M 2 2 V(*) + 4V(t). (2.29) 



2.4 Spherical fields 

We now examine the results of the previous section in the context of spherical 
field theory. We start with the spherical partial wave expansion, 

<\> = E(=0,l,-X n =O r .(E m =-„,.. Ii ^,n,m(i)ii,n, m (^ <f), (2.30) 

where Y^ m ^ n are four- dimensional spherical harmonics satisfying 

/ d 3 Q y / * n / >m /(6 ) , if), ip)Yi jn>m (9, if), If) = SvjSni^Sm^m, (2.31) 

10 If our theory contained bare masses similar arguments would apply for the 

2 

infrared limit /i 2 , mf —>■ 0, with — j- fixed. 

11 There may be additional logarithmic factors but this does not matter for our 
purposes here. 
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1 l,n,mK u i 



Vw) = (-l) m nn,-m(0,^)- (2-32) 

The explicit form of Y^ m ^ n can be found in . The integral of the free massless 
Lagrange density in terms of spherical fields is 



fd*tC = f~dt { £ [(-l) m <f>l,n, 
V l,m,n 



d_£_d_ t 

dt 2 dt 



-U(l + 2) 



'Ln,m 



(2.33) 



It can be shown that the process of angle smearing the field </>(t) is equivalent to 
multiplying 4>i^ m (t) by an extra factor s z M (t) where 

_ 2Mt[(t+2)sin(^)-t B in(^)] . . 

S l W- ia+l)G+2)[2-Mt S in(^)j • 

For large /, sf 4 ^) diminishes as Z~ 2 , and so the correlator receives an extra suppression 
of Z~ 4 . We will later use this result to estimate the contribution of high spin partial 
waves. The regularization of our correlator can be implemented in our Lagrange 
density by dividing factors of s^it), 



■ > Ln,—m 



JLi-JL _ 1/(7 + 21 
at. 2 at 2H l t 



(2.35) 



[( S f(t))-V*,n,-J 



at 2 at 



-li(i+2) [(sfwrvw] 



We now include the interaction and counterterms. We first define 



l3,n 3 ,m 3 ;l4, 714,7714 

(0, ^ ^)^4 ,714,7714 



(2.36) 



We can write the full functional integral as 

| X?0exp [/d 4 f£] ex I (n^ m ^m) ex P Uo°° dt (^i + L 2 + 2*)] , (2.37) 
where 

Li = E [(^(OrVL-J [|tI - 1^ + 2)] [(^W)-V^nj] , (2-38) 

Z, 771,71 



12 



TA ] deserves credit as the first discussion of radial (or covariant Euclidean) quan- 



tization, an important part of the spherical field formalism. 
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l,m,n 



'-11 



— C- 



(0) 



V«iV</> c V(/>V</> 



ti,(0) 



+c4 iv^(/+2)+t3(M : 



a. £ 

at 2 

2-(0) 



_a 



Ln.m 



(2.39) 



,(oh 



E[ii,mi,ni: 
\_l 3 ,m 3 ,n 3 : 



l2,m 2 ,n 2 
£4, 7714,714 



ll,mi,niTl2 ,m2 ,712 ^l3,m 3 ,713 $14,1714,114 " 



(2.40) 



We have used primes in preparation for redefining the fields, 



(^(*))~VLm = <f>l,n,m. 



(2.41) 



The Jacobian of this transformation is a constant (although infinite) and can be 
absorbed into the normalization of the functional integral. Now the Lagrangian Li 
has the usual free-field form in terms of </>z, n ,m while L2 and L3 are now functions of 

,n,m- 

With M serving as our ultraviolet regulator, the contribution of high-spin partial 
waves decouples for sufficiently large spin I. We can estimate the order of magni- 
tude of this contribution in the following manner. We first identify t~ 1 l (where t 
is the characteristic radius we are considering) as an estimate of the magnitude of 
the tangential momentum, px- For px 3> M t^ 1 our correlator scales as By 
dimensional analysis, a diagram with Nl loops and Nj internal lines will receive a 
contribution from partial waves with spin > I of order 



M 4 



\4JV L 



M 4 



(rH) 



(2.42) 



2.5 One- loop examples 

We will devote the remainder of our discussion to computing one-loop spherical 
Feynman diagrams as a check of our results. Our calculations are done both numer- 
ically and analytically. The diagrams we will consider are shown in Figures ^]2| and 
We start with the two-point function in Figure [2T2"|. The amplitude can be 
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written as t 3 B(t) where 

Constants of proportionality are not important here and so we will define B(t) to 
be equal to the right side of ( |2.43|) . Our results tell us that if we choose our mass 
counterterms appropriately, the combination 

B(t)+M^ + ^ (2.44) 

should be independent of t, or more succinctly, 

B(t) + £cg (2.45) 

is independent of t. Let us first check this analytically. In the absence of a high-spin 
cutoff, we can explicitly calculate the sum in (|2.43| ): 



where 



In the limit M — > oo, 



B(t) = ± + b(t) (2.46) 



B{t) - £ + |M 2 . (2.48) 



We conclude that c^ 2 = —1 and Bit) + ^c^/ is in fact translationally invariant. 
In Figure |2.4| we have plotted B{t) — p-, computed numerically for various values 



of the high-spin cutoff J max - We have also plotted the limiting values b{t) and |M 2 . 
In our plot t is measured in units of m~ l and B(t) — h is in units of m 2 , where m is 
an arbitrary mass scale such that M = 3m. As expected, the errors are of size ^jp-. 
There is clearly a deviation from |M 2 for t < M _1 but the integral of the deviation 
is negligible as M — ► oo. 

We now turn to the four-point function in Figure |2.3| . The amplitude can be 
written as t^C^ti, £2) where 

C(h,h) cx E.nJ-^Brp^ \^d{t 2 - tl ) + - 1 2 )} 2 . (2.49) 



u 2 t 1 
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Again constants of proportionally are not important and so we will define C(t x , t 2 ) to 
be equal to the right side of ( |2.49| ). We can write C(ti,t 2 ) in terms of the regulated 
correlator G^fti, fr^P^l 

C(t u t 2 ) oc J d%d 3 i 2 [G M 2(r i ,f 2 )] 2 oc J d% [G M ^{ti, t 2 )] 2 . (2.50) 

Since the coupling constant counterterm 

$U\t x -t 2 ) (2.51) 

is translationally invariant, the amplitude by itself should be translationally invariant. 
Let us define 

r d%e~^-^ [G M2 (t x ,t 2 )] 2 = f(f 



so that 



d%e^ [G M 2(ti,t 2 )] = e^f(f). 
Integrating over ii, we find 

dt 2 t 2 2 J x (pt 2 )C(t x ,t 2 ) oc ±J x (pt x )f{f). 



Let us define 



We now check that in fact 



C(t) = fdt 2 t 2 2 J x (pt 2 )C(t,t 2 ). 



C(t) oc ±J x (pt x ). 

In the absence of a high-spin cutoff, we find that C(t) is given byQ 
C{t) = ±J x {pU 



where the ellipsis represents terms which vanish as M — > 00 and 

-1/2 



c = 324 



ii I (sin k— k cos k) 4 1 \ . / ,JZ~( s i 

aK \ 4F3 324fe j + y i 2 a/C_ 



sin fc— fc cos fc) 4 
iF3 



(2.52) 

(2.53) 

(2.54) 

(2.55) 
(2.56) 

(2.57) 
(2.58) 



In Figure ^75] we plot C{t) for different values of the high-spin cutoff J max as well as 
13 We recall that the regulated correlator goes with i>\ nm rather than (f)i >n>m . But 

this is not important here since 4>' Q00 = 0o,o,o- 

14 This calculation is somewhat lengthy. Details can be obtained upon request from 

the authors. 



Chapter 2: Renormalization in spherical field theory 



20 




the large-M limit value 



Ci(f) = ±Ji(ptil 



+ C 



(2.59) 



In our plot t is measured in units of p 1 and M = 3p. From (|2.42j) the expected error 



is of size ^fp-. We see that the data is consistent with the results expected. Again 



the deviation for t < M 1 integrates to a negligible contribution as M — > oo. 



2.6 Summary 

We have examined several important features of non-perturbative renormaliza- 
tion in the spherical field formalism and answered the three questions posed in the 
introduction. Ultraviolet divergences can be cancelled by a finite number of local 
counterterms in a manner such that the renormalized theory is translationally in- 
variant. Using angle-smearing regularization we find that the counterterms for 4 
theory in four dimensions can be parameterized by five unknown constants as shown 
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in ( 2.29Q . Aside from our remarks about Ward identity constraints in gauge theories, 
the extension to other field theories is straightforward. We hope that these results 
will be useful for future studies of general renormalizable theories by spherical field 
techniques. 



Chapter 3 

The massless Thirring model in 
spherical field theory i 



3.1 Introduction 



The massless Thirring model |58j is an exactly soluble system of a single self- 
interacting massless fermion in two dimensions. There are a number of solutions 
in the literature based on properties of the Euler-Lagrange equations and fermion 
currents or bosonization techniques f 



5|,|g|ig|BE3Sa»S»Ei@- Given 



its simplicity and solubility, the model has become a popular testing ground for new 
ideas and methods in field theory. 

From a computational point of view, however, the massless Thirring model still 
presents a significant challenge. In the lattice field formalism, the difficulties are 
due to the appearance of fermion doubler states and singular inversion problems 
associated with integrating out massless fermions. In this work we use the model to 
illustrate new non-perturbative methods in the spherical field formalism [35, |36], [5], [38 



The techniques we present are quite general and can also be applied to other modal 
expansion methods such as periodic field theory 



As noted in [36], we will not need to worry about fermion doubling. This is 



X N. Salwen, D. Lee, Phys. Lett. B468 (1999) 118. 
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true for any modal field theory and follows from the fact that space is not discretized 
but retained as a continuous variable. Since our model is not super-renormalizable 
we will use a procedure called angle-smearing, a regularization method designed for 
spherical field theory |3^| . With angle-smearing regularization and a small set of 
local counterterms, we are able to remove all ultraviolet divergences in a manner such 
that the renormalized theory is finite and translationally invariant. Comparison of 
our results with the known Thirring model solution will serve as a consistency check 
for our regularization and renormalization procedures. 

The organization of this paper is as follows. We begin with a short summary 



of the massless Thirring model, following the solution of Hagen [20|. Using angle- 
smearing regularization we obtain the spherical field Hamiltonian and construct a 
matrix representation of the Hamiltonian. We reduce the space of states using a two- 
parameter auxiliary cutoff procedure. In this reduced space we compute the time 
evolution of quantum states using a fourth-order Runge-Kutta-Fehlberg algorithm. 
We calculate the two point correlator for several values of the coupling and find 
agreement with the known analytic solution. 



3.2 The model 

We start with a list of our notational conventions. Our analysis will be in two- 
dimensional Euclidean space, and we use both cartesian and polar coordinates, 



t = (t cos 9, t sin 6) = (x,y). 



The components of the spinors ip and ijj are written as 



(3.1) 



■0 



(3.2) 



Our representation for the Dirac matrices is 



7 = iff, 



(3.3) 
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and so 7 satisfies 

{ 7 \y} = -2^, ij = 1,2. (3.4) 
The massless Thirring model is formally defined by the Lagrange density 

£ = ^7-W-U-J, (3.5) 



where j is the fermion vector current. Johnson |3(J emphasized the importance of 
defining the regularized current precisely, and this was further clarified by the work of 



Sommerfield |55| and Hagen p0| . We will use a regularization technique, introduced 
3"8|, called angle-smearing. We define the regularized current as 



111 



3 = I tyal^s - Trtfip s i/j s ]) , (3.6) 

where 



M^ 9 ) = / deip(t,6 + e). (3.7) 

We identify the radial variable t as our time parameter, and our definition of the 
current is local with respect to t. 

Hagen [2(J described the solution of the Thirring model in the Hamiltonian for- 
malism with currents defined as products of the canonical operators at equal times. 
Though our equal-time surface is curved, the curvature of the integration segment in 
Q3.7D scales as 4? while the ultraviolet divergences in this model are only logarithmic 
in M. In the M — > 00 limit we therefore recover the standard results. As discussed in 
20|] , there exists a one parameter class of solutions to the Thirring model depending 



on the preferred definition of the regularized vector and axial vector currents. We 



will use the conventions used in |30| and |pq| , which in Hagen's notation corresponds 
with the parameter values £ = r\ — ^. With this choice the Hamiltonian density 
takes the form 

H = H free + ^(t • j) 2 + ^ c (§- j) 2 , (3.8) 

where 

c=± (3.9) 

The only counterterms needed in this model are wavefunction renormalization 
counterterms, a result of our careful definition for the regularized interaction in 
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As in [^UJ we calculate correlation functions using an unrenormalized Hamiltonian. 
The divergent wavefunction normalizations will appear simply as overall factors in 
the correlators. 



3.3 Spherical field Hamiltonian 



In this section we derive the form of the spherical field Hamiltonian. We first 
expand the fermion current in terms of components of the spinors, 

e~ ie 



9 ■ ips^ips = 4> s 



3 t0 

e~* 
-e ie 



^ = e~ l W s - e*W s . 



The anti-commutators of the regulated fields areQ 



= \ {W f^Le^de = A(t)e ie 



(3.10) 



(3.11) 



(3.12) 



where 



= \ ) 2 /5e- i W (fe = A(t)e 



From the anti-commutation relations, the t component of the current is 

t-J = i [ie- ie (^ l s - #1) + ie*$W, - 
= ie~ ie ^ l s + ie*^jV>I - iA(t), 

and so 

{t ■ j) 2 = A(t) [e-^ l s + e i6 ^jl] - 2#Vi^J - A 2 (t). 
Similarly we find 

■ j) 2 = A(t) [e~ i6 ^i + e»$ij>l] - 2<MVW 



(3.13) 
(3.14) 

(3.15) 



(3.16) 



(3.17) 



'Our definition of the Euclidean fermion fields and anti-commutation relations 



follows the conventions of 14 . 
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The Hamiltonian can now be written as 



H = H free + Jd6t{^ [A(t) [e-*&tyi + e w ^ s } - 2^i</M] } . (3.18) 



We have omitted the constant term proportional to A 2 (t). 
Let us define the partial wave modes 

1>n(t) = ^kf d6e- me ^(i), <Mt) = ^J d6e^ s (t), 
It is straightforward to show that for n ^ 0, 

AM = s iWMt) ^ >B (t) = 7Si(t). 

\Mt) \ Mt J 



Sin (l?t) - 



\Mt) 

At this point it is convenient to rescale ■?/>, 



ft = w n - 



In terms of the partial waves, 



^A(t)sin(^)sin(S±i) » t T / / J. 
\Mt)\ Mt ) 



-IE 



n — 



6 ^(*) sin (^) sin (W)\^_^T 



(3.19) 
(3.20) 

(3.21) 



We can extend this result to the case n = using the convenient shorthand 

= 1 . (3.22) 



(3.23) 



(3.24) 



E 



,/,T sm(^) smt^- ) sm(^- ) sm(^) 
n\ Vn 2 +1 V-n 3 -l rn 4 ? "i \7 "2+1 V7 "3+ 1 YTJhA 

\MtJ\ Mt J\ Mt )\Mt) 



where 



2c 
1-c 2 



(3.25) 



Since b is the parameter appearing in the Hamiltonian, it is somewhat more convenient 
to express c in terms of b, 

c= yJ±E^ (326) 
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Let us define the ladder operators^ 



a -n> a -n — Vru 1p- n -l 
1 ^ -ri + l,^n. 



a n+li a n+l 



(3.27) 
(3.28) 



These operators satisfy the anti-commutation relations 



i a i a it\ = ; a T a Tt\ = < 5 



1 12 ■ 



(3.29) 



with all other anti-commutators equal to zero. We can now recast the Hamiltonian 

as 



H 



IE 



n + 



^tA(t) S m(^)sm(^) 



\ Mt ) \ Mt J 



E 

-ni+n2+n3— n4=0 



i u ni u, n2 U, n3 U, n 4 



fan 4 + «n a l) 

L r )sm(-^-)sm(- 

Aft A Mt A Mt J \ Mt J 



(3.30) 



sin ( ^ir- ) sin ( wi ) sin ( I Mt i ) sin ( jh ) 



We will implement a high spin cutoff by removing terms in the interaction contain- 
ing operators a n , a^, a n , or a n ' for |n| > J ma x- This has the effect of removing high 
spin modes, which correspond with large tangential momentum states. We should 
emphasize, however, that J max is an auxiliary cutoff. It does not play a role in the 
regularization scheme since the interactions have already been rendered finite using 
angle-smearing. The contribution of these high spin modes is negligible so long as 



(3.31) 



where t is the characteristic radius of the process being measured. Returning back to 
( |3.121 ) and ( 3.1 3| ) and removing the contribution of these partial waves, we find that 
A{t) is replaced by 



A., 



n=-J u 



\Mt) V Mt ) 



(3.32) 



follows: ai,a^ — n ^ + ' " ' - " " ] ~ r 



This notation is slightly different from that used in [36 . The translation is as 
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Let \0)j ree be the ground state of the free massless fermion Hamiltonian.^ For 



n > 0, we find 



4 \°)free = 4 W free = ^ (3-33) 

and for n < 0, 

oj? |0) /ree = oJJ |0) /ree = 0. (3.34) 
It is convenient to define the normal-ordered products 

i + i J a^al for n > T+ t J a)^a) n for n > 

■ a "<- = \ . : « : = < t t . (3-35) 

[ — a^a^' for n < [ — a^a^ for n < 0. 

The ordering for other operators is immaterial since the anti-commutators are zero. 
We can now rewrite H in terms of normal-ordered products, 



H=(^ + 0(J-l)) (a^ai + a^ai) 



(3.36) 



-ni+n 2 +n 3 — n 4 =0 



t ■ u n 1 u 'n2 U 'nz u 'n i - 



sm (^) sm (irn:) sm ( mt) sm (i77) 

( n l- L \( n 2\( n 3~ 1 \( 

V Mt )\ Mt ) \ Mt ) \ Mt ) 



There is an 0(J~^ X ) term due to a small asymmetry in our cutoff procedure with 
respect to the two boundaries — J max and J max -Q We will neglect this term in the 
limit J max — > oo. 



3.4 Two-point correlator 

We wish to study the properties of the two-point correlator. The massless Thirring 
model is invariant under the discrete transformation 



fiAjHt) - -^ T (*~), (3.37) 

4 The ground state of the free massless Hamiltonian is actually degenerate due to 
s-wave excitations, but this is remedied by taking the m — > limit of the massive 
fermion theory. 

5 If desired we could eliminate this term and the asymmetry by a slight change in 
the angle-smearing procedure for ■0. 
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as well as the transformation 



From these we deduce 



(0| T [^(t)^(0)} |0> = <0| T $H*V(0)] |0) = 



and 



(0\T[^{t, 9)^(0)] |0) = (0\T[^(t, -9)^(0)] |0). 
It therefore suffices to consider just the correlator on the left side of ( 3.40Q . 



In the limit M — > oo the form of the correlator is given by 

(0\T[^(t, 9)^(0)] |0> = ^{k{c)M)^t^~'~ 



(3.38) 



(3.39) 



(3.40) 



(3.41) 



where k(c) is a dimensionless parameter. Standard analytic methods do not yield 
a simple closed form expression for k(c). We will therefore extract k(c) from the 
computed value of the correlator at a specific renormalization point t = io-Q 
We define 

f(t) = (0\T[i>\(t)4(0)] |0>. (3.42) 

Since 



<0| T [^(t, 9)^(0)} |0> = £ <0| T te(f)^(0)l |0) , 



we conclude that 



-2c^ -1-c^ 



f(t) = (k(c)M)^?t^ r . 



(3.43) 



(3.44) 



We now compute f(t) using the spherical field Hamiltonian. We first need a ma- 
trix representation for the Grassmann ladder operators. We will use tensor products 
of the 2x2 identity matrix and Pauli matrices: 

6 In some regularization schemes k(c) can be calculated analytically p3| |[44|| , and 
it may be worthwhile to use these techniques in future work. In this first analysis, 
however, we prefer to present a more straightforward and typical example of the 
angle-smearing regularization method. 
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a l n = °z + 1, (3.45) 

fc— Jxu&Xf ^max ^ — ^max ,TlH~l %—Tl 1 , </ max 

4 = 0- Z ® fax + fry) 10 1- 

^ — ^maxj^"l"l % — Tl 1 , t/max % = J max i t/max 

The representations for a$ and are defined by the conjugate transposes of these 
matrices. We can now calculate the correlator f(t) using the relation [[36 



/(«) = lim lim ^{-S^^ff-M-I^^^ (3 . 46) 

A straightforward calculation of ( |3.46| ), however, is rather inefficient. There are sev- 
eral techniques which we will first use to simplify the calculation. 

The time evolution of the system at large t is dominated by the contribution of 
the ground state or, more precisely, the adiabatic flow of the t-dependent ground 
state. As discussed in |35|, |36|, |[ a similar phenomenon occurs at small t, due to the 
divergence of energy levels near t = 0. It is therefore not necessary to compute the 
full matrix trace in the numerator and denominator of ( 3.46|) . It is instead sufficient 



to compute the corresponding ratio for a single matrix element. After making this 
reduction, we can then go a step further and eliminate states which do not contribute 
to the matrix element. 

The high spin parameter J max was used to remove high spin modes with \n\ > J max - 
This, however, is not a uniform cutoff in the space of states and most of the remaining 
states are still high kinetic energy states. Although none of the individual modes are 
energetic, many of the modes can be simultaneously excited. Let us define and 

to be bit switches, 1 or 0, depending on whether or not the corresponding mode is 
excited. Let us also define a cutoff parameter K max . We will remove all high kinetic 
energy states such that 

En [N ( N i + N D] > ^max- (3.47) 

For consistency K max should be about the same size as J max - 
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3.5 Results 



The CPU time and memory requirement for calculating (|3.46|) scales linearly with 
the number of transitions in H (i.e., non-zero elements in our matrix representation). 
In Table 1 we have shown the number of states and transitions for different values of 



Table 1 





2 


4 


6 


8 


10 


12 


states 


6 


40 


210 


920 


3600 


13000 


transitions 


38 


500 


4200 


26000 


1.4E5 


3.9E5 



We have calculated f(t) for J max < 12 and several values of the coupling b. The total 
run time was about 100 hours on a 350 MHz PC with 256 MB RAM. 

The matrix time evolution equations in ( |3.46| ) were computed using a fourth-order 
Runge-Kutta-Fehlberg algorithm. We have set 



(3.48) 



We will use the notation /j max (t) to identify the corresponding result for a given value 
of Jmax- In Figure |3J] we have plotted fj mea .(t) for b = 1 and J max = 4,8, 12. We 
have scaled t and f(t) in dimensional units chosen such that M = 3. For finite J max 



we expect deviations from the J max — > oo limit to be of size 0{y^-). The curves 



shown in Figure |3.1| appear consistent with this rate of convergence. 

We can extrapolate to the limit J max — > oo using the asymptotic form 



(3.49) 



For b = 0,0.5,1.0,3.0 and M = 3 we have calculated f(t) using this extrapolation 
technique for J max = 10 and 12.[] The results are shown in Figure 3j . For comparison 
we have plotted the analytic solution 



f A (t) = (k(c)M)T^ t - 



(3.50) 



7 Both our results and the analytic solution are even in b, and so we consider only 
positive values. 
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The relation between b and c can be found in (|3.25|) and (|3.26|) . The parameter 
k(c) is fitted to the value of the correlator at the renormalization point t = 0.6.0 The 
agreement appears quite good. Some deviations from the analytic solution are due 
residual terms, which were left out of the derivation of (|3.50|). These 



to 01 



i 



APt 2 > 



effects are significant in the small t region, t < M l , especially for larger values of b. 
The values we find for k(c) are shown in Table 2.0 



Table 2 



b 


0.5 


1.0 


3.0 


k(c) 


1.77 


1.77 


1.68 



We can compare our results at small b with a simple perturbative calculation. Evaluat- 
ing the corresponding regulated two-loop diagram we obtain, for small b, k(c) ~ 1.75. 

8 The relative error is expected to be small in the vicinity of this point. 

9 In some regularization schemes k(c) can be shown to be independent of the cou- 
pling. Our regularization method seems to be rather close to this, with only a slow 
variation with respect to the coupling strength. 
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This appears consistent with the results in Table 2. 

3.6 Summary 

We derived the angle-smeared spherical Hamiltonian for the massless Thirring 
model and constructed an explicit matrix representation. We discarded negligible 
high energy states using auxiliary cutoff parameters J max and K max . In this reduced 
space we computed the time evolution of quantum states and calculated the two-point 
correlator for several values of the coupling. The results of our computation are in 
close agreement with the known analytic solution. In addition to demonstrating 
new computational methods, our analysis also serves as a consistency check of the 
regularization and renormalization methods introduced in |3^| . 

We believe that this work represents a significant new direction in the non- 
perturbative computation of fermion dynamics. Future work will study the ap- 
plication of these methods to systems of interacting bosons and fermions. 



Chapter 4 



Modal expansions and 
non-perturbative quantum field 
theory in Minkowski spaced 

4.1 Introduction 

Modal expansion methods have recently been used to study non-perturbative phe- 
nomena in quantum field theory [ |35| , |36| , [| [38|]. Modal field theory, the name for the 
general procedure, consists of two main parts. The first is to approximate field the- 
ory as a finite-dimensional quantum mechanical system. The second is to analyze the 
properties of the reduced system using one of several computational techniques. The 
quantum mechanical approximation is generated by decomposing field configurations 
into free wave modes. This technique has been explored using both spherical partial 
waves (spherical field theory fl35|, |36| , ||, |5l] ) and periodic box modes (periodic field 
theory @). 

Having reduced field theory to a more tractable quantum mechanical system, we 
have several different ways to proceed. Boson interactions in Euclidean space, for 
example, can be modeled using the method of diffusion Monte Carlo. In many situ- 

X N. Salwen, D. Lee, Phys. Rev. D62 (2000) 025006. 
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ations, however, Monte Carlo techniques are inadequate. These include unquenched 
fermion systems, processes in Minkowski space, and the phenomenology of multi- 
particle states. Difficulties arise when the functional integral measure cannot be 
treated as a probability distribution or when information must be extracted from 
excited states obscured by dominant lower lying states. Fortunately there are sev- 
eral alternative methods in the modal field formalism which avoid these problems. 
Matrix Runge-Kutta techniques were introduced in [ET]] as a method for calculating 
unquenched fermion interactions. Here we discuss a different approach, one which di- 
rectly calculates the spectrum and eigenstates of the Hamiltonian. For this approach 
it is essential that the Hamiltonian is time-independent, and so we will consider peri- 
odic rather than spherical field theory. As we demonstrate, these methods naturally 
accommodate the study of multi-particle states and Minkowskian dynamics. 

We apply the spectral approach to 1 + 1 dimensional 4 theory in a periodic box 
and calculate the real and imaginary parts of the <fi propagator. Some interesting 
properties of <fi\ theory such as the phase transition at large coupling were already dis- 
cussed within the modal field formalism using Euclidean Monte Carlo techniques 
The purpose of this analysis is of a more general and exploratory nature. Our aim is 
to test the viability of modal diagonalization techniques for quantum field Hamilto- 
nians. We would like to know whether we can clearly see multi-particle phenomena, 
the size of the errors and computational limitations with current computer resources, 
and how such methods might be extended to more complicated higher dimensional 
field theories. 

The spectral method presented in the first part of our analysis is similar to the 
work of Brooks and Frautschi |27], ^6| ,Q who considered a 1 + 1 dimensional Yukawa 
model in a periodic box and deserves credit for the first application of diagonalization 
techniques using plane wave modes in a periodic box. Our calculations are also 



similar in spirit to diagonalization-based Hamiltonian lattice formulations |23| , |24fl and 
Tamm-Dancoff light-cone and discrete light-cone quantization |17], [18, |2|, || . There 



are, however, some differences which we should mention. As in [27, 26 



we are 



2 We thank the referee of the original manuscript for providing information on this 
reference. 
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using a momentum lattice rather than a spatial lattice. We find this convenient to 
separate out invariant subspaces according to total momentum quantum numbers. 
Since we are using an equal time formulation our eigenvectors are not boost invariant 
as they would be on the light cone. Also we are using a simple momentum cutoff 
scheme rather than a regularization scheme which includes Tamm-Dancoff Fock-space 
truncation. As a result our renormalization procedure is relatively straightforward, 
but we will have to confront the problem of diagonalizing large Fock spaces from the 
very beginning. In the latter part of the paper we mention current work on quasi- 
sparse eigenvector methods which can handle even exceptionally large Fock spaces. 
Despite the differences among the various diagonalization approaches to field theory, 
the issues and problems discussed in our analysis are of a general nature. We hope 
that the ideas presented here will be of use for the various different approaches. 



4.2 Spectral method 

The field configuration <ft in 1 + 1 dimensions is subject to periodic boundary 
conditions <p(t, x — L) = <ft(t, x + L). Expanding in terms of periodic box modes, we 
have 

0(t, x) = v /T J- Mt)e mnx/L . (4.1) 

n=0,±l,... 

The sum over momentum modes is regulated by choosing some large positive number 
Amax and throwing out all high-momentum modes <p n such that \n\ > iV max . In this 
theory renormalization can be implemented by normal ordering the 4 interaction 
term. After a straightforward calculation (details are given in |[43||), we find that the 
counterterm Hamiltonian has the form 

Ji E ( 4 - 2 ) 

ti = iV max ,A^ max 

where 

b= E 3b ^ = V ^ + ^ ( 4 - 3 ) 

TL= -/Vniax , -/V m ax 
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We represent the canonical conjugate pair <p n and using the Schrodinger opera- 
tors q n and Then the functional integral for <p A theory is equivalent to that 
for a quantum mechanical system with Hamiltonian 



H = 



£ [-Iafci + l(^-t)?-n?n] (4.4) 



~^ 4!2L ^ ^ (Zm 9«2 (JVi3 ■ 



17>i — ^max^max 

ni+n2+ri3+n4=0 

We now consider the Hilbert space of our quantum mechanical system. Given d, 
an array of non-negative integers, 

d = {d- Nmax , ■ ■ ■ d NmBX } , (4.5) 

we denote Pd(q) as the following monomial with total degree \d\, 

n=-A r max,A'max n 

We define to be a Gaussian of the formP] 



G((q) = ] [ exp 



9 



-ngnA/C 2 +" 2 ^ 2 /^ 2 



(4.7) 



£ is an adjustable parameter which we will set later. Any square- integrable function 
ip(q) can be written as a superposition 

= 5>p d ( 9 )Gr f (g). (4.8) 

d 

In this analysis we consider only the zero-momentum subspace. We impose this 
constraint by restricting the sum in ( |4.8j ) to monomials satisfying 



J2nd n = 0. (4.9) 



3 G((q) has been defined such that G^(q) is the ground state of the free theory. 



Chapter 4-' Modal expansions 



39 



We will restrict the space of functions if)(q) further by removing high energy states 
in the following manner. Let 

k(d) = \n\ d n . (4.10) 



k(d) was first introduced in |51| and provides an estimate of the kinetic energy associ- 
ated with a given state. Let us define two auxiliary cutoff parameters, K max and D max . 
We restrict the sum in ( |4.8| ) to monomials such that k(d) < K max and \d\ < D max . 
We will refer to the corresponding subspace as Vx max ,D max - The cutoff K max removes 
states with high kinetic energy and the cutoff D max eliminates states with a large 
number of excited modes.0 We should stress that K max and D max are only auxil- 
iary cutoffs. We increase these parameters until the physical results appear close to 
the asymptotic limit K max , D max — > oo. In our scheme ultraviolet regularization is 
provided only by the momentum cutoff parameter N max . 

Our plan is to analyze the spectrum and eigenstates of H restricted to this approx- 
imate low energy subspace, Vftr max ,D max - For any fixed L and iV max , H is the Hamil- 
tonian for a finite-dimensional quantum mechanical system and the results should 
converge in the limit K max , D max — > oo. We obtain the desired field theory result by 
then taking the limit L, jV '" ax — > oo. 



4.3 Results 

We have calculated the matrix elements of H restricted to V 7 ft- maxi D max using a 
symbolic differentiation-integration algorithm^ and diagonalized H, obtaining both 
eigenvalues and eigenstates. Let A be the full propagator, 

A(p 2 ) = i [ d 2 x e lv » xV (0| T [0(^)0(0)] |0) . (4.11) 



4 In the case of spontaneous symmetry breaking, the broken symmetry of the vac- 
uum may require retaining a large number of go modes. This however is remedied 
by shifting the variable, q' Q = qo — (0| go |0). 

5 A11 codes can be obtained by request from the authors. 
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We have computed A by inserting our complete set of eigenstates (complete in 
^max.-Dmax)- Let A mp be the multi-particle contribution to A, 

A mp (p 2 ) = A(p 2 )-A pole (p 2 ), (4.12) 

where A po i e is the single-particle pole contribution. We are primarily interested in 
A mp , a quantity that cannot be obtained for p 2 > using Monte Carlo methods. 
Since the imaginary part of A po i e is a delta function, it is easy to distinguish the 
single-particle and multi-particle contributions in a plot of the imaginary part of A. 
The real part of A, however, is dominated by the one-particle pole. For this reason 
we have chosen to plot the real part of A mp rather than that of A. 

Although we have referred to multi-particle states, it should be noted that in 
our finite periodic box there are no true continuum multi-particle states. Instead 
we find densely packed discrete spectra with level separation of size ~ L^ 1 which 
become continuum states in the limit L — > oo. We can approximate the contribution 
of these L — > oo continuum states by a simple smoothing process. We included a 
small width r ~ L^ 1 to each of the would-be continuum states and averaged over a 
range of values for L. For the results we report here we have averaged over values 
L = 2.07T, 2.l7r, • • -2.871". For convenience all units have been scaled such that fi = 1. 

The parameter ( was adjusted to reduce the errors due to the finite cutoff values 
-f^max and -Dmax- Since the spectrum of H is bounded below, errors due to finite K milx 
and -D max generally drive the estimated eigenvalues higher. One strategy, therefore, 
is to optimize ( by minimizing the trace of H restricted to the subspace V^ maxi £> max . 
The approach used here is a slight variation of this — we have minimized the trace of 
H restricted to a smaller subspace Vk< d< C Vr- „ n „. The aim is to accelerate 

1 max ' max max ; J -' max 

the convergence of the lowest energy states rather than the entire space V^ maxi £> max . 
Throughout our analysis we used K' max , D' max = 8,3. 



For = 0.50 we have plotted the imaginary part of A in Figure |4.1| and the real 
part of A mp in Figure [4.2| . The value 4 = 0.50 is above the threshold for reliable 
perturbative approximation^ (|y < 0.25) but below the critical value at which — > —<p 
symmetry breaks spontaneously (4 ~ 2.5). The contribution of the one-particle 

6 For momenta \p 2 \ > 1. 
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Im [A (p 2 




Figure 4.1: Imaginary part of A(jo 2 ) for 4 = 0.50 and several values for 



^maxj ^maxi D miix 



Re[A mp (p 2 ) ] V4! = 0.50 




Figure 4.2: Real part of A mp (p 2 ) for ^ = 0.50 and several values for iV max , K m£kX , D 
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state appears near p 2 = (0.93) 2 and the three-particle threshold is at approximately 
p 2 = (2.9) 2 . We have chosen several different values for iV max , K max , -D max to show 
the convergence as these parameters become large. The plot for iV max , K max , D max = 
9, 19, 7 appears relatively close to the asymptotic limit. The somewhat bumpy texture 
of the curves is due to the finite size of our periodic box and diminishes with increasing 
L. From dimensional power counting, we expect errors for finite iV max to scale as 
iV~ ax . Assuming that K m£ix and -D max also function as uniform energy cutoffs, we 
expect a similar error dependence - and it appears plausible from the results in Figures 



[4. 1| and fO|. A more systematic analysis of the errors and extrapolation methods for 
finite iV max , -Kmax, -D max , and L, will be discussed in future work. 



In Figures [4.3| and [4.4| we have compared our spectral calculations with the two- 



Im[A(p^) ] 

4xlCT 6 



V4! = 0.01 

N K D =9197 



10 



spectral calculation 
two-loop result 
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Figure 4.3: Imaginary part of A(jo 2 ) for ^ = 0.01 and comparison with the two-loop 
result 



loop perturbative result for 4 = 0.01. We have used iV max , K max , D max = 9, 19, 7, 
and the agreement appears good. For small 4 the propagator has a very prominent 
logarithmic cusp at the three-particle threshold, which can be seen clearly in Figures 
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A./ 4! = 0.01 
Re [A m „ (p )] m k D =9197 




spectral calculation 
two-loop result 



Figure 4.4: Real part of A mp (p 2 ) for 4 = 0.01 and comparison with the two-loop 
result 
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^and O 



In Figures fO>] and £0] we have compared results for 4 = 0.25, 0.50, 1.00. We 




Figure 4.5: Imaginary part of A(p 2 ) for ^ = 0.25,0.50, 1.00 



have again used iV max , i^ max , -D max = 9, 19, 7. In contrast with the quadratic scaling 
in the perturbative regime, the results here scale approximately linearly with 4. An 
interesting and perhaps related observation is that the magnitude of the multi-particle 
contribution to A remains small (< 0.003) even for the rather large coupling value 



4.4 Limitations and new ideas 

We now address the computational limits of the diagonalization techniques pre- 
sented in this work. These techniques are rather straightforward and can in principle 
be generalized to any field theory. In practise however the Fock space Vx max ,D max 
becomes prohibitively large, especially for higher dimensional theories. The data in 
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Re [A mp (p 2 )] N max ,K max ,D„ 
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Figure 4.6: Real part of A mp (p 2 ) for | = 0.25, 0.50, 1.00 



Figures [O] and |4.2j and crosschecks with Euclidean Monte Carlo results[] suggest that 
for iV max = 9 and L = 2.07T, • • • 2.87T our spectral results with K max , D mauX = 19, 7 and 
4 < 1 are within 20% of the K max , -D max — > oo limit. In this case Vft- maXi £) max is a 2036 
dimensional space and requires about 100 MB of RAM using general (dense) matrix 
methods. 

Sparse matrix techniques such as the Lanczos or Arnoldi schemes allow us to 
push the dimension of the Fock space to about 10 5 states. This may be sufficient 
to do accurate calculations near the critical point 4 ~ 2.5 for larger values of L 
and A^ max . It is, however, near the upper limit of what is possible using current 
computer technology and existing algorithms. Unfortunately field theories in 2 + 1 
and 3 + 1 dimensions will require much larger Fock spaces, probably at least 10 12 and 
10 18 states respectively. In order to tackle these larger Fock spaces it is necessary 
to venture beyond standard diagonalization approaches. The problem of large Fock 
spaces (3>10 6 states) is beyond the intended scope of this analysis. But since it 



7 See ESI for a discussion of these methods. 
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is of central importance to the diagonalization approach to field theory we would 
like to briefly comment on current work being done which may resolve many of the 
difficulties. The new approach takes advantage of the sparsity of the Fock-space 
Hamiltonian and the approximate (quasi-)sparsity of the eigenvectors. A detailed 
description will be provided in a future publication [j39| . 

We start with some observations about the eigenvectors of the <Pi +1 Hamiltonian 
for iV max = 9, L = 2.57T and K max , D max = 19, 7. To make certain that we are probing 
non-perturbative physics we will set ^ = 2.5, the approximate critical point value. 
We label the normalized eigenvectors as |t>o), • • , ascending in order with respect 
to energy. We also define \bo), \bi),- ■ ■ as the normalized eigenvectors of the free, 
non-interacting theory. For any v,i we know 

EK fe >*>i 2 = L ( 4 - 13 ) 

3 

Let us define |||ui)|L n ) as the partial sum 

llk)ll ( n)= E KWI 2 , (4.i4) 

k=l,---n 

where the inner products have been sorted from largest to smallest 

\(b h \vi)\ > \(b h \vi)\ >■•• . (4.15) 
Table 1 shows Hl^i)!!^ for several eigenvectors and different values of n. 

Table 1 





IK)ll(n) 


n= 10 


n = 20 


n = 40 


n = 80 




^o) 


0.75 


0.84 


0.90 


0.94 






0.89 


0.92 


0.95 


0.97 




V5) 


0.87 


0.91 


0.94 


0.96 




Vio) 


0.77 


0.86 


0.90 


0.94 



Despite the non-perturbative coupling and complex phenomena associated with the 
phase transition, we see from the table that each of the eigenvectors can be approxi- 
mated by just a small number of its largest Fock-space components. We recall that 
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the Fock space for this system has 2036 dimensions. The eigenvectors are therefore 
quasi-sparse in this space, a consequence of the sparsity of the Hamiltonian. If we 
write the Hamiltonian as a matrix in the free Fock-space basis, a typical row or col- 
umn contains only about 200 non-zero entries, a number we refer to as A^ ransition . In 
|39| we show that a typical eigenvector will be dominated by the largest V A^ transition 
elements.^ The key point is that V A^ transition is quite manageable — on the order of 
10 3 and 10 5 for 2 + 1 and 3 + 1 dimensional field theories respectively. Although 
the size of the Fock space for these systems are enormous, the extreme sparsity of 
the Hamiltonian suggests that the eigenvectors can be approximated using current 
computational resources. 

With this simplification, the task is to find the important basis states for a given 
eigenvector. Since the important basis states for one eigenvector are generally differ- 
ent from that of another, each eigenvector is determined independently. This provides 
a starting point for parallelization, and many eigenvectors can be determined at the 
same time using massively parallel computers. In |5{J we present a simple stochas- 
tic algorithm where the exact eigenvectors act as stable fixed points of the update 
process. 



4.5 Summary 

We have introduced a spectral approach to periodic field theory and used it to 
calculate the propagator in 1 + 1 dimensional 4 theory. We find that the straightfor- 
ward application of these methods with existing computer technology can be useful for 
describing the multi-particle properties of the theory, information difficult to obtain 
using Euclidean Monte Carlo methods. However the extension to higher dimensional 
theories is made difficult by the large size of the corresponding Fock space. As a pos- 
sible solution to this problem, we note that each eigenvector of the <fi\ + i Hamiltonian 
can be well-approximated using relatively few components and discuss some current 



There are some special exceptions to this rule and they are discussed in [39]. But 
these are typically not relevant for the lower energy eigenstates of a quantum field 
Hamiltonian. 
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work on quasi-sparse eigenvector methods. 



Chapter 5 

Quasi sparse methods in the 
diagonalization of quantum field 



5.1 Introduction 

Most computational work in non-perturbative quantum field theory and many 
body phenomena rely on one of two general techniques, Monte Carlo or diagonaliza- 
tion. These methods are nearly opposite in their strengths and weaknesses. Monte 
Carlo requires relatively little storage, can be performed using parallel processors, and 
in some cases the computational effort scales reasonably with system size. But it has 
great difficulty for systems with sign or phase oscillations and provides only indirect 
information on wavefunctions and excited states. In contrast diagonalization meth- 
ods do not suffer from fermion sign problems, can handle complex-valued actions, 
and can extract details of the spectrum and eigenstate wavefunctions. However the 
main problem with diagonalization is that the required memory and CPU time scales 
exponentially with the size of the system. 

In view of the complementary nature of the two methods, we consider the combi- 




^ean Lee, N. Salwen, Dan Lee, Phys. Lett. B 503(2001) 223-235 
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nation of both diagonalization and Monte Carlo within a computational scheme. We 
propose a new approach which takes advantage of the strengths of the two compu- 
tational methods in their respective domains. The first half of the method involves 
finding and diagonalizing the Hamiltonian restricted to an optimal subspace. This 
subspace is designed to include the most important basis vectors of the lowest energy 
eigenstates. Once the most important basis vectors are found and their interactions 
treated exactly, Monte Carlo is used to sample the contribution of the remaining basis 
vectors. By this two-step procedure much of the sign problem is negated by treating 
the interactions of the most important basis states exactly, while storage and CPU 
problems are resolved by stochastically sampling the collective effect of the remaining 
states. 

In our approach diagonalization is used as the starting point of the Monte Carlo 
calculation. Therefore the two methods should not only be efficient but work well 
together. On the diagonalization side there are several existing methods using Tamm- 
Dancoff truncation | i?fl , similarity transformations |T7|, 19], density matrix renormal- 
ization group |BC, [J]]], or variational algorithms such as stochastic diagonalization 



25| , fL0|. However we find that each of these methods is either not sufficiently gen- 



eral, not able to search an infinite or large dimensional Hilbert space, not efficient 
at finding important basis vectors, or not compatible with the subsequent Monte 
Carlo part of the calculation. The Monte Carlo part of our diagonalization/Monte 
Carlo scheme is discussed separately in a companion paper jHJ. In this paper we 
consider the diagonalization part of the scheme. We introduce a new diagonaliza- 
tion method called quasi-sparse eigenvector (QSE) diagonalization. It is a general 
algorithm which can operate using any basis, either orthogonal or non-orthogonal, 
and any sparse Hamiltonian, either real, complex, Hermitian, non- Hermit ian, finite- 
dimensional, or infinite-dimensional. It is able to find the most important basis 
states of several low energy eigenvectors simultaneously, including those with identi- 
cal quantum numbers, from a random start with no prior knowledge about the form 
of the eigenvectors. 

Our discussion is organized as follows. We first define the notion of quasi-sparsity 
in eigenvectors and introduce the quasi-sparse eigenvector method. We discuss when 
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the low energy eigenvectors are likely to be quasi-sparse and make an analogy with An- 
derson localization. We then consider three examples which test the performance of 
the algorithm. In the first example we find the lowest energy eigenstates for a random 
sparse real symmetric matrix. In the second example we find the lowest eigenstates 
sorted according to the real part of the eigenvalue for a random sparse complex non- 
Hermitian matrix. In the last example we consider the case of an infinite-dimensional 
Hamiltonian defined by 1 + 1 dimensional 4 theory in a periodic box. We conclude 
with a summary and some comments on the role of quasi-sparse eigenvector diago- 
nalization within the context of the new diagonalization/Monte Carlo approach. 

5.2 Quasi-sparse eigenvector method 

Let \ei) denote a complete set of basis vectors. For a given energy eigenstate 



we define the important basis states of \v) to be those |ej) such that for fixed normal- 
izations of \v) and the basis states, |q| exceeds a prescribed threshold value. If \v) 
can be well-approximated by the contribution from only its important basis states we 
refer to the eigenvector \v) as quasi-sparse with respect to |ej). 

Standard sparse matrix algorithms such as the Lanczos or Arnoldi methods allow 
one to find the extreme eigenvalues and eigenvectors of a sparse matrix efficiently, 
without having to store or manipulate large non-sparse matrices. However in quan- 
tum field theory or many body theory one considers very large or infinite dimensional 
spaces where even storing the components of a general vector is impossible. For these 
more difficult problems the strategy is to approximate the low energy eigenvectors 
of the large space by diagonalizing smaller subspaces. If one has sufficient intuition 
about the low energy eigenstates it may be possible to find a useful truncation of the 
full vector space to an appropriate smaller subspace. In most cases, however, not 
enough is known a priori about the low energy eigenvectors. The dilemma is that 
to find the low energy eigenstates one must truncate the vector space, but in order 
to truncate the space something must be known about the low energy states. 




(5.1) 
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Our solution to this puzzle is to find the low energy eigenstates and the appro- 
priate subspace truncation at the same time by a recursive process. We call the 
method quasi-sparse eigenvector (QSE) diagonalization, and we describe the steps 
of the algorithm as follows. The starting point is any complete basis for which the 
Hamiltonian matrix Hij is sparse. The basis vectors may be non-orthogonal and/or 
the Hamiltonian matrix may be non- Hermit ian. The following steps are now iterated: 

1. Select a subset of basis vectors {e^ , • • • , ej n } and call the corresponding subspace 
S. 

2. Diagonalize H restricted to S and find one eigenvector v. 

3. Sort the basis components of v according to their magnitude and remove the 
least important basis vectors. 

4. Replace the discarded basis vectors by new basis vectors. These are selected 
at random according to some weighting function from a pool of candidate basis 
vectors which are connected to the old basis vectors through non-vanishing 
matrix elements of H. 

5. Redefine S as the subspace spanned by the updated set of basis vectors and 
repeat steps 2 through 5. 

If the subset of basis vectors is sufficiently large, the exact low energy eigenvectors 
will be stable fixed points of the QSE update process. We can show this as follows. 
Let \i) be the eigenvectors of the submatrix of H restricted to the subspace S, where 
S is the span of the subset of basis vectors after step 3 of the QSE algorithm. Let 
\Aj) be the remaining basis vectors in the full space not contained in S. We can 
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represent H as 



Ai 






A 2 



(l|#|Aa) 
(2|^|A) 



(1\H\A 2 ) 
(2\H\A 2 ) 



{A X \H\\) {A X \H\2) 
(A 2 \H\1) (A 2 \H\2) 



E-X Al (Ax\H\A 2 ) 
(A 2 \H\A 1 ) E-X M 



(5.2) 



We have used Dirac's bra-ket notation to represent the terms of the Hamiltonian 
matrix. In cases where the basis is non-orthogonal and/or the Hamiltonian is non- 
Hermitian, the meaning of this notation may not be clear. When writing (A± \ H\l), 
for example, we mean the result of the dual vector to \Ai) acting upon the vector 
H |1). In (|5.2p we have written the diagonal terms for the basis vectors \Aj) with 
an explicit factor E. We let |1) be the approximate eigenvector of interest and have 
shifted the diagonal entries so that Ai = 0. Our starting hypothesis is that |1) is close 
to some exact eigenvector of H which we denote as |lf u ii). More precisely we assume 
that the components of |lf u n) outside S are small enough so that we can expand in 
inverse powers of the introduced parameter E. 
We now expand the eigenvector as 

1 



1 1 full) 



d 2 E~ l + ■■■ 

d A E^ + ... 
c' A2 E-' + ... 



(5.3) 



and the corresponding eigenvalue as 

Ami = ^[E 1 + ■ ■ 



In ( |5.3| ) we have chosen the normalization of |lf u ii) such that (1 |lf u ii) 
eigenvalue equation 

H | lfuLl) = Af u u |lfull) 



(5.4) 
1. From the 

(5.5) 
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we find at lowest order 

We see that at lowest order the component of |lf u ii) in the \Aj) direction is independent 
of the other vectors \Ajr). If |1) is sufficiently close to |lf u ii) then the limitation that 
only a fixed number of new basis vectors is added in step 4 of the QSE algorithm is 
not relevant. At lowest order in E~ l the comparison of basis components in step 3 
(in the next iteration) is the same as if we had included all remaining vectors \Aj) 
at once. Therefore at each update only the truly largest components are kept and 
the algorithm converges to some optimal approximation of |lf u ii). This is consistent 
with the actual performance of the algorithm as we will see in some examples later. 
In those examples we also demonstrate that the QSE algorithm is able to find several 
low energy eigenvectors simultaneously. The only change is that when diagonalizing 
the subspace S we find more than one eigenvector and apply steps 3 and 4 of the 
algorithm to each of the eigenvectors. 



5.3 Quasi-sparsity and Anderson localization 

As the name indicates the accuracy of the quasi-sparse eigenvector method de- 
pends on the quasi-sparsity of the low energy eigenstates in the chosen basis. If 
the eigenvectors are quasi-sparse then the QSE method provides an efficient way to 
find the important basis vectors. In the context of our diagonalization/Monte Carlo 
approach, this means that diagonalization does most of the work and only a small 
amount of correction is needed. This correction is found by Monte Carlo sampling 
the remaining basis vectors, a technique called stochastic error correction pD| . If 



however the eigenvectors are not quasi-sparse then one must rely more heavily on the 
Monte Carlo portion of the calculation. 

The fastest and most reliable way we know to determine whether the low energy 
eigenstates of a Hamiltonian are quasi-sparse with respect to a chosen basis is to use 
the QSE algorithm and look at the results of the successive iterations. But it is 
also useful to consider the question more intuitively, and so we consider the following 
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example. 

Let if be a sparse Hermitian 2000 x 2000 matrix defined by 



H jk = log(j) • 6 jk + x jk ■ M. 



(5.7) 



where j and k run from 1 to 2000, xjf. is a Gaussian random real variable centered 
at zero with standard deviation x rms = 0.25, and Mjk is a sparse symmetric matrix 
consisting of random 0's and l's such that the density of l's is 5%. The reason 
for introducing the log(j) term in the diagonal is to produce a large variation in 
the density of states. With this choice the density of states increases exponentially 
with energy. Our test matrix is small enough that all eigenvectors can be found 
without difficulty We will consider the distribution of basis components for the 



eigenvectors of H. In Figure 5.1 we show the square of the basis components for 




200 400 600 800 1000 1200 1400 1600 1800 2000 

sorted basis vector |e> 



Figure 5.1: Distribution of basis components for an eigenvector where the spacing 
between consecutive levels is AE = 0.13x rTTls . 



a given low energy eigenvector \v) . The basis components are sorted in order of 
descending importance. The ratio of AE, the average spacing between neighboring 
energy levels, to x rms is 0.13. We see that the eigenvector is dominated by a few 
of its most important basis components. In Figure we show the same plot for 
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Figure 5.2: Distribution of basis components for an eigenvector where AE 
0.04Lr rms . 



another eigenstate but one where the spacing between levels is three times smaller, 
AE/x rms = 0.041. This eigenvector is not nearly as quasi-sparse. The effect is even 
stronger in Figure |5^, where we show an eigenvector such that the spacing between 
levels is AE/x rms = 0.024. 

Our observations show a strong effect of the density of states on the quasi-sparsity 
of the eigenvectors. States with a smaller spacing between neighboring levels tend 
to have basis components that extend throughout the entire space, while states with 
a larger spacing tend to be quasi-sparse. The relationship between extended versus 
localized eigenstates and the density of states has been studied in the context of 



Anderson localization and metal-insulator transitions [y, 56, |63|, 28]. The simplest 
example is the tight-binding model for a single electron on a one-dimensional lattice 
with Z sites, 

H £ (l > •/) 0' ■ £ 0/ ./) (/ • (5-8) 
\j) denotes the atomic orbital state at site j, dj is the on-site potential, and tjf is 
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Figure 5.3: Distribution of basis components for an eigenvector where AE = 
0.024x rms . 

the hopping term between nearest neighbor sites j and j'. If both terms are uniform 
(dj = d, tjji = t) then the eigenvalues and eigenvectors of H are 



where n = 1, • ■ ■ , Z labels the eigenvectors. In the absence of diagonal and off- 
diagonal disorder, the eigenstates of H extend throughout the entire lattice. The 
eigenvalues are also approximately degenerate, all lying within an interval of size 
At. However, if diagonal and/or off-diagonal disorder is introduced, the eigenvalue 
spectrum becomes less degenerate. If the disorder is sufficiently large, the eigenstates 
become localized to only a few neighboring lattice sites giving rise to a transition of 
the material from metal to insulator. 

We can regard a sparse quantum Hamiltonian as a similar type of system, one 
with both diagonal and general off-diagonal disorder. If the disorder is sufficient 
such that the eigenvalues become non-degenerate, then the eigenvectors will be quasi- 




(5.9) 




(5.10) 
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sparse. We reiterate that the most reliable way to determine if the low energy states 
are quasi-sparse is to use the QSE algorithm. Intuitively, though, we expect the 
eigenstates to be quasi-sparse with respect to a chosen basis if the spacing between 
energy levels is not too small compared with the size of the off-diagonal entries of the 
Hamiltonian matrix. 



5.4 Finite matrix examples 

As a first test of the QSE method, we will find the lowest four energy states of the 
random symmetric matrix H defined in (|5.7|) . So that there is no misunderstanding, 
we should repeat that diagonalizing a 2000 x 2000 matrix is not difficult. The purpose 
of this test is to analyze the performance of the method in a controlled environment. 
One interesting twist is that the algorithm uses only small pieces of the matrix and 
operates under the assumption that the space may be infinite dimensional. A sample 
MATLAB program similar to the one used here has been printed out as a tutorial 
example in |R1 . 



The program starts from a random configuration, 70 basis states for each of the 
four eigenvectors. With each iteration we select 10 replacement basis states for each of 



the eigenvectors. In Figure 5.4 we show the exact energies and the results of the QSE 



method as functions of iteration number. In Figure [5.5| we show the inner products 
of the normalized QSE eigenvectors with the normalized exact eigenvectors. We note 
that all of the eigenvectors were found after about 15 iterations and remained stable 
throughout successive iterations. Errors are at the 5 to 10% level, which is about the 
theoretical limit one can achieve using this number of basis states. The QSE method 
has little difficulty finding several low lying eigenvectors simultaneously because it 
uses the distribution of basis components for each of the eigenvectors to determine 
the update process. This provides a performance advantage over variational-based 
techniques such as stochastic diagonalization in finding eigenstates other than the 
ground state. 

As a second test we consider a sparse non-Hermitian matrix with complex eigen- 
values. This type of matrix is not amenable to variational-based methods. We will 
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Figure 5.4: Comparison of the four lowest exact energies Ei and QSE results E® SE 
as functions of iteration number. 



Eigenvector inner products 
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Figure 5.5: Inner products between the normalized exact eigenvectors |i>j) and the 



QSE results 



as functions of iteration number. 
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find the four eigenstates corresponding with eigenvalues with the lowest real part for 
the random complex non-Hermitian matrix 

H' jk = (l + i-c jk )H jk . (5.11) 

Hjk is the same matrix used previously and uniform random variable dis- 

tributed between —1 and 1. As before the program is started from a random con- 
figuration, 70 basis states for each of the four eigenvectors. For each iteration 10 



replacement basis vectors are selected for each of the eigenvectors. In Figure |5.6| the 

Complex eigenvalues 



80 - 




3 - 1 - 5 lm(E) 
Re(E) 



Figure 5.6: Comparison of the four lowest exact eigenvalues (sorted by real part) 
and QSE results E® SE in the complex plane as functions of iteration number. 



exact eigenvalues and the results of the QSE run are shown in the complex plane as 
functions of iteration number. In Figure |5.7| we show the inner products of the QSE 
eigenvectors with the exact eigenvectors. All of the eigenvectors were found after 
about 20 iterations and remained stable throughout successive iterations. Errors 
were again at about the 5 to 10% level. 
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Figure 5.7: Inner products between the normalized exact eigenvectors and the 



QSE results 
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as functions of iteration number. 



5.5 



theory in 1 + 1 dimensions 



We now apply the QSE method to an infinite dimensional quantum Hamiltonian. 
We consider 4 theory in 1 + 1 dimensions, a system that is familiar to us from previous 
studies using Monte Carlo ||3 and explicit diagonalization [52|. The Hamiltonian 
density for 4 theory in 1 + 1 dimensions has the form 

W=Hff)' + M2)' + '^ + 

where the normal ordering is with respect to the mass \i. We consider the system in 
a periodic box of length 2L. We then expand in momentum modes and reinterpret 



the problem as an equivalent Schrodinger equation ||43|| . The resulting Hamiltonian 
is 



n n 



(5.12) 



+ 



4!2L J j QniQ^QnaQn* 

ni+n2+n3+ri4=0 
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where 



UJ n (fi) 



and b(fj,) is the coefficient for the mass counterterm 



(5.13) 



(5.14) 



It is convenient to split the Hamiltonian into free and interacting parts with respect 
to an arbitrary mass //: 



n n 

H = H free + ±£ (/ - f - It 1 ) 

n 

rti+n2+ri3+n4=0 



(5.15) 



(5.16) 



// is used to define the basis states of our Fock space. Since H is independent of //, 
we perform calculations for different // to obtain a reasonable estimate of the error. 
It is also useful to find the range of values for // which maximizes the quasi-sparsity 
of the eigenvectors and therefore improves the accuracy of the calculation. For 
the calculations presented here, we set the length of the box to size L = 57r/i _1 . We 
restrict our attention to momentum modes q n such that \n\ < iV max , where iV max = 20. 
This corresponds with a momentum cutoff scale of A = 4/i. 

To implement the QSE algorithm on this infinite dimensional Hilbert space, we 
first define ladder operators with respect to //, 



4 (//) = . 1 

The Hamiltonian can now be rewritten as 



q n u n {n') + ^ 



(5.17) 
(5.18) 



H = J2^>la n + \{f - ^ - ^py 



~ 192L / j 

ni+n2+"3+«4=0 



("ni+nlnj {an 2 +al n2 ) (an 3 +at ra3 ) (ar^+a^J 



(5.19) 



Chapter 5: Quasi sparse methods 



63 



In ( 5.19|) we have omitted constants contributing only to the vacuum energy. We rep- 



resent any momentum-space Fock state as a string of occupation numbers, |o_7v max , • • • , 0jv max )> 
where 

a n a n |o_7v max , ■ • • , o Nin J = o n \o_ Nmax , ■ ■ • , o Nrn J . (5.20) 

From the usual ladder operator relations, it is straightforward to calculate the matrix 
element of H between two arbitrary Fock states. 

Aside from calculating matrix elements, the only other fundamental operation 
needed for the QSE algorithm is the generation of new basis vectors. The new states 
should be connected to some old basis vector through non-vanishing matrix elements 
of H. Let us refer to the old basis vector as |e). For this example there are two 
types of terms in our interaction Hamiltonian, a quartic interaction 

^ ] ^ a rii + (^ a ri2 + a -n2^j ^ a "3 a -n^j (^ a ~n 1 ~n 2 -n 3 + 0-n 1 + r! , 2 + r! ,3^ , (5.21) 

and a quadratic interaction 

( a ~n + 4.) («». + a-n) ■ (5-22) 

n 

To produce a new vector from |e) we simply choose one of the possible operator 
monomials 

^rt2 a n ^a— ni _ n 2 — r«3 j o_ ni (3 n2 ct n3 ct_ ni _ n2 _ n3 , , (5.23) 
ct_ n (3 n , a n a_ n , 

and act on |e). Our experience is that the interactions involving the small momentum 
modes are generally more important than those for the large momentum modes, a 
signal that the ultraviolet divergences have been properly renormalized. For this 
reason it is best to arrange the selection probabilities such that the smaller values of 
| rail, | ri2 1, |ras| and |ra| are chosen more often. 

For each QSE iteration, 50 new basis vectors were selected for each eigenstate and 
250 basis vectors were retained. The results for the lowest energy eigenvalues are 



shown in Figure [5.8|. The error bars were estimated by repeating the calculation for 
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Figure 5.8: Energy eigenvalues of as functions of the coupling constant. 



different values of the auxiliary mass parameter //. 

From prior Monte Carlo calculations we know that the theory has a phase transi- 



tion at 



2.5/i 2 corresponding with spontaneous breaking of the 



—6 reflection 



symmetry. In the broken phase there are two degenerate ground states and we refer 



to these as the even and odd vacuum states. In Figure |5.8| we see signs of a second 
order phase transition near 4 « 2.5/i 2 . Since we are working in a finite volume 
the spectrum is discrete, and we can track the energy eigenvalues as functions of the 
coupling. Crossing the phase boundary, we see that the vacuum in the symmetric 
phase becomes the even vacuum in the broken phase while the one-particle state in 
the symmetric phase becomes the odd vacuum. The energy difference between the 
states is also in agreement with a Monte Carlo calculation of the same quantities. 
The state marking the two-particle threshold in the symmetric phase becomes the 
one-particle state above the odd vacuum, while the state at the three-particle thresh- 
old becomes the one-particle state above the even vacuum. These one-particle states 
should be degenerate in the infinite volume limit. One rather unusual feature is the 
behavior of the first two-particle state above threshold in the symmetric phase. In 
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the symmetric phase this state lies close to the two-particle threshold. But as we 
cross the phase boundary the state which was the two-particle threshold is changed 
into a one-particle state. Thus our two-particle state is pushed up even further to 
become a two-particle state above the even vacuum and we see a pronounced level 
crossing. 

We note that while the one-particle mass vanishes near the critical point, the 
energies of the two-particle and three-particle thresholds reach a minimum but do 
not come as close to zero energy. It is known that this model is repulsive in the 
two-particle scattering channel. In a large but finite volume the ground state and 
one-particle states do not feel significant finite volume effects. The two-particle state 
at threshold, however, requires that the two asymptotic particles be widely separated. 
In our periodic box of length 2L the maximal separation distance is L and we expect 
an increase in energy with respect to twice the one-particle mass of size ~ V(L), 
where V is the potential energy between particles. Likewise a three-particle state 
will increase in energy an amount ~ 3V(2L/3). Our results indicate that finite 
volume effects for the excited states are significant for this value of L. 

5.6 Summary 

We have proposed a new approach which combines both diagonalization and 
Monte Carlo within a computational scheme. The motivation for our approach 
is to take advantage of the strengths of the two computational methods in their re- 
spective domains. We remedy sign and phase oscillation problems by handling the 
interactions of the most important basis states exactly using diagonalization, and 
we deal with storage and CPU problems by stochastically sampling the contribution 
of the remaining states. We discussed the diagonalization part of the method in 
this paper. The goal of diagonalization within our scheme is to find the most im- 
portant basis vectors of the low energy eigenstates and treat the interactions among 
them exactly. We have introduced a new diagonalization method called quasi-sparse 
eigenvector diagonalization which achieves this goal efficiently and can operate using 
any basis, either orthogonal or non-orthogonal, and any sparse Hamiltonian, either 
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real, complex, Hermitian, non-Hermitian, finite-dimensional, or infinite-dimensional. 
Quasi-sparse eigenvector diagonalization is the only method we know which can ad- 
dress all of these problems. 

We considered three examples which tested the performance of the algorithm. 
We found the lowest energy eigenstates for a random sparse real symmetric matrix, 
the lowest eigenstates (sorted according to the real part of the eigenvalue) for a 
random sparse complex non-Hermitian matrix, and the lowest energy eigenstates 
for an infinite-dimensional Hamiltonian defined by 1 + 1 dimensional <fr theory in a 
periodic box. 

We regard QSE diagonalization as only a starting point for the Monte Carlo 
part of the calculation. Once the most important basis vectors are found and their 
interactions treated exactly, a technique called stochastic error correction is used to 
sample the contribution of the remaining basis vectors. This method is introduced 
inlOj. 

Acknowledgments We thank P. van Baal, H. S. Seung, H. Sompolinsky, and M. 
Windoloski for useful discussions. 



Chapter 6 

Introduction to stochastic error 
correction methods^ 



6.1 Introduction 



In [[39] a new approach was proposed for finding the low-energy eigenstates of very 
large or infinite-dimensional quantum Hamiltonians. This proposal combines both 
diagonalization and Monte Carlo methods, each being used to solve a portion of the 
problem for which the technique is most efficient. The first part of the proposal is to 
diagonalize the Hamiltonian restricted to a subspace containing the most important 
basis vectors for each low energy eigenstate. This may be accomplished either through 
variational techniques or an ab initio method such as quasi-sparse eigenvector (QSE) 
diagonalization. The second step is to include the contribution of the remaining basis 
vectors by Monte Carlo sampling. The use of diagonalization allows one to consider 
systems with fermion sign oscillations and extract information about wavefunctions 
and excited states. The use of Monte Carlo provides tools to handle the exponential 
increase in the number of basis states for large volume systems. 



The first half of this proposal was discussed in [39|. An adaptive diagonalization 



algorithm known as the quasi-sparse eigenvector method was introduced to find the 



1 D. Lee, N. Salwen, M. Windoloski, Phys. Lett. B 502(2001) 329-337 
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most important basis vectors for each low energy eigenstate. This technique is espe- 
cially valuable when little is known about the low energy states. It is also the only 
method available which can handle non-orthogonal bases, non-Hermitian Hamiltoni- 
ans, infinite dimensional systems, and which can find several low energy states with 
like quantum numbers simultaneously. In this paper we discuss the second half of 
the diagonalization/Monte Carlo scheme. We introduce several new Monte Carlo 
techniques which we call stochastic error correction (SEC). There are two general 
varieties of stochastic error correction, methods based on a series expansion and those 
which are not. The series method starts with an eigenvector of the Hamiltonian re- 
stricted to some starting subspace and then includes the contribution of the remaining 
basis states as terms in an ordered expansion. The idea is to form a perturbative 
expansion centered around a good non-perturbative starting point. 

As an example of a non-series method we discuss a technique called the stochastic 
Lanczos method. This method again starts with eigenvectors of a Hamiltonian sub- 
matrix. Using these as starting vectors, we define Krylov vectors, \j) ,H \j) , H 2 \j) ■ ■ ■ , 
similar to standard Lanczos diagonalization. The new ingredient is that matrix ele- 
ments between Krylov vectors, \ H n \ j) , are computed using matrix diffusion Monte 
Carlo. Since the method does not rely on a series expansion, it has the advantage 
that the starting vectors need not be close to the exact eigenvectors. 

One can generate a large class of stochastic error correction methods based on 
other non-series algorithms, various ways of resumming the series expansion, or com- 
binations of the two techniques. In this introductory paper we concentrate on describ- 
ing the basic principles and implementation of the series and non-series approaches. 
We also present three test problems which demonstrate the potential of the new ap- 
proach for a range of different problems. In the first example we determine the low 
energy spectrum of 02+i theory using QSE diagonalization and first order corrections 
using the series method. In the second example we find the low energy spectrum 
of compact U(l) in 2 + 1 dimensions using the stochastic Lanczos method. In the 
last example we find the ground state of the 2 + 1 dimensional Hubbard model using 
QSE diagonalization and first order series stochastic error correction. In each case 
we compare with published results in the literature. We conclude with a summary 
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and some general comments on the new computational scheme. 



6.2 Series method 



Let \i) be the eigenvectors of a Hamiltonian H restricted to some subspace S. 
Let \Aj) be the remaining basis vectors in the full space not contained in S. We can 
represent H as 



Ai 






A 2 



{A x \H\l) (A 1 \H\2) 
(A 2 \H\1) {A 2 \H\2) 



(l\H\A x ) ll\H\A 2 ) 

(2\H\A X ) (2\H\A 2 ) 

E-X Al (A X \H\A 2 ) 
{A 2 \H\A X ) E-\ A2 



(6.1) 



We have used Dirac's bra-ket notation to represent the terms of the matrix. In cases 
where the basis is non-orthogonal or the Hamiltonian is non-Hermitian, the precise 
meaning of terms such as (A x \ H |1) is the action of the dual vector to \A X ) upon the 
vector H |1). We have written the diagonal terms for the basis vectors \Aj) with an 
explicit factor E for reasons to be explained shortly. 

Let us assume that |1) is close to some exact eigenvector of H which we denote as 
| lfuii) . More concretely we assume that the components of |lf u ii) outside S are small 
enough so that we can expand in inverse powers of the introduced parameter E. In 
order to simplify the expansion we choose to shift the diagonal entries so that Ai = 0. 

The series method of stochastic error correction is based on the E~ x expansion, 



I lfuii) oc 



c'E- 1 + c"E- 2 + 



c> Al E^ + c" Al E-* + 
c'mE' 1 + c" A E-* + 



(6.2) 
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Afiui = a;^- 1 + \';e- 2 ■■■ . (6.3) 

It is convenient to choose the normalization of the eigenvector such that the |1) com- 
ponent remains 1. The convergence of the expansion is controlled by the proximity 
of |1) to |lf u n). If |1) is not at all close to |lf u ii) then it will be necessary to use 
a non-series method such as the stochastic Lanczos method discussed in the next 
section. 

At first order in E' 1 we find 



c A . = 



(Aj\H\l) 
A A, 



(6.4) 



(l|H|^-)(A,-|H|r) 
Aa, 



c '. = jl m\A k )(A k \H\i) 

3 A, \ A 



At second order the contributions are 



1^1 k 



(6.5) 
(6.6) 



c „ = (A 3 \H\A k )(A k \H\l) ST^Y^ {A 3 \H\l)(l\H\A k )(A k \H\l) , g ^ 



\" = \^ (MH\Aj)(AmAk)(A k \H\l) _Spspsp (llHlA^jA^H^illHlA^jAklHll) , g g , 

1 / J / J Aa^Aaj. / j / j / j Aa^AjAaj, 

j fe^j j Z^l fe 



EE 



(m|H|^>(Aj|g|^ fc >(A fc |g|l> 



+EEE 



AmAA J A(AA fc 



(6.9) 



E 



Aa, 



E (m|H|A fc )(A fc |g|l) 
AA t 



These contributions can be calculated by straightforward Monte Carlo sampling. All 
that is required is an efficient way of generating random basis vectors \ A k ) with known 
probability rates. Let P(-A tr iai) denote the probability of selecting | Atrial) on a given 
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trial. If for example we are calculating the first order correction to the eigenvalue, 
then we have 

K = _Y^ m\AMA m i) (61Q) 

j 3 

lim 1 ^triaKi)) (^trial(i) pi 1 ) 

i=l,-,N 

6.3 Stochastic Lanczos 

We now consider a method called stochastic Lanczos which does not require the 
starting vectors to be close to exact eigenvectors of H. This is essential if the eigen- 
vectors of H are not quasi-sparse and require extremely large numbers of basis states 
to represent accurately. 

Let V be the full Hilbert space for our system. As in the previous section let S be 
the subspace over which we have diagonalized H exactly. Let P$ be the projection 
operator for S and let Xj and \j) be the eigenvalues and eigenvectors of H restricted 
to S so that 

P s HP s \j) = X j \j). (6.11) 

Let Z be an auxiliary subspace, one which contains S but excludes very high-energy 
states. Let Pz be the projection operator for Z. We will choose Z such that 
P Z HP Z is bounded above. Let a be a real constant which is greater than the 
midpoint of the minimum and maximum eigenvalues of PzHPz- As n — > oo the 
operator [Pz(H — a)Pz] n maps any given state in Z to the corresponding lowest- 
energy eigenvector of PzHPz with non-zero overlap. 

The stochastic Lanczos method uses the operators [Pz{H — a)P z ] n to approximate 
the low-energy eigenvalues and eigenvectors of PzHP z . The goal is to diagonalize 
if in a subspace spanned by vectors 

\d,j) = [P z (H-a)P z ] d \j), (6.12) 

for several values of d and j. This requires calculating (d',j'\ d,j) and (d',j'\ H \d,j). 
If our Hamiltonian matrix is Hermitian, both of these terms can be written in the 
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general form 

(j'\[P z (H-a)P z ] n \j). (6.13) 
Therefore it suffices to determine the matrix 

A n = P s [P Z (H - a)P z ] n P s . (6.14) 

For non-orthogonal bases and non-Hermitian Hamiltonians, the only change is that 
we use vectors 

[P z (H-a)P z ] d \j) (6.15) 
to generate approximate right eigenvectors of H and vectors in the dual space 

(j\[P z (H-a)P z ] d (6.16) 

to produce approximate left eigenvectors. Adding and subtracting Ps(H — a)Ps, we 
can rewrite 

A n = P s [Pz(H - a)P z - P S (H - a)P s + P S (H - a)P s ] n P s . (6.17) 
A n can now be evaluated recursively as 

A n+1 = B n+1 + B ™( H ~ a )A n - m , (6.18) 

m=0,--- ,n 

where 

B n = P s [Pz(H - a)P z - P S (H - a)P s ] n P s . (6.19) 

The components of B n are computed by matrix diffusion Monte Carlo. One 
could also directly evaluate the components of A n . However the calculation for B n 
eliminates the need to sample the matrix Ps(H — a)Ps, which is already known. Any 
general matrix product M^M^ ■ ■ ■ is a sum of degree n monomials, 

[M«M< 2 >...M%= M % M £l-- M tlk- (6-20) 



We can interpret ( |6.20| ) as a sum over paths through the set of basis vectors of Z 1 

\j) _, \ k ) _, . . . _, _> \ k ) , (6.21) 
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with an associated weight M^M-^ ■ ■ ■ M^ lk . The components of B n are sampled 



using ensembles of random walkers. We refer the interested reader to |32| for a review 
of methods in diffusion Monte Carlo. 

We end the section with a discussion of the fermion sign problem. The sign 
problem is a general issue for any Monte Carlo calculation. For a system with sign 
oscillations the evaluation of a Euclidean-time Green's function involves sums 

( 6 - 22 ) 



with the property that 



^ iXi - exp(-c-y-T), (6.23) 



where V is the volume, T is the Euclidean time, and c is a positive constant. We 
will refer to this term as the cancellation ratio. The exponential dependence on V 
and T makes computations difficult even for small systems. 

The sign problem will affect the calculation of B n in the stochastic Lanczos algo- 
rithm and terms in the series method discussed in the previous section. The effect 
however is different from the sign problem in typical Monte Carlo Green's function 
calculations. Stochastic error correction is a calculation of eigenvalues and eigenvec- 
tors rather than a sampling of the partition function or the time evolution of a given 
initial state. Therefore the quantity of interest is not exp(-HT) but the action of 
H or H n on approximate eigenvectors of H. Due to homogeneity in H the explicit 
volume dependence does not appear in the cancellation ratio. Instead we find 

j^', ~exp(-fc-w), (6.24) 

Li N 

where k is a positive constant. The sign problem will return if the starting point of 
the SEC calculation is very poor and it becomes necessary to use n such that k ■ n 
is large. However in many cases k ■ n can be kept small even for large n since the 
most important part of the Hamiltonian, PsHPs, is diagonalized exactly. In short 
the sign problem is less severe because stochastic error correction uses the result of 
subspace diagonalization as its starting point. 
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6.4 4 theory in 2 + 1 dimensions 

The first example we consider is 4 theory in 2 + 1 dimensions near the 0^—0 
symmetry restoration phase transition. We will use QSE diagonalization and the 
series version of stochastic error correction to probe the low energy spectrum of the 
theory on both sides of the phase transition. 

In 1 Magruder demonstrated the existence of a phase transition in 02+i by 



extending Chang's duality argument for <fi\ + i. The statement of the main result is 
as follows. Consider the two Lagrange densities 

£+ = \d v <j>dr<j> - \^ - ^ + \8l +( j> 2 (6.25) 
C_ = Id^d'cj) + ±/x 2 _0 2 - + |^_0 2 . (6.26) 

The counterterm 5 2 + is defined so that in the C + system the self-energy graphs 
vanish at zero-momentum up to two-loop order. By shifting the field 



<P = <P , + \j'^Y ( 6 - 27 ) 
we note that the same counterterm S 2 (same mass dependence but /x + replaced by 
//_) is also sufficient to renormalize By equating C + and £_ we obtain a duality 
constraint between the two theories. One feature of this constraint is that the g —>■ oo 
limit of £_ is mapped to the g — > limit of C + . Therefore whose reflection 
symmetry — > —<p is broken at small g, must eventually reach the symmetric phase 
for sufficiently large coupling.^ 

The £_ phase transition was studied using quasi-sparse eigenvector diagonaliza- 
tion with stochastic error correction. Quantities such as the critical coupling, critical 
exponents, and the low lying energy spectrum were studied and, where possible, com- 
pared with Monte Carlo results. A full discussion methods and results are presented 
in p2"| . We will very briefly summarize some of the results below. 

The two spatial dimensions of our system are taken to be a periodic box of size 

2L by 2L. We will use the modal field formalism to describe the Hamiltonian for the 

2 The g — ► oo limit of C + is mapped to the g — > oo limit of £_ and so there is no 
analogous argument for a phase transition in C + . Numerical calculations indicate 
that there is no phase transition for C + . 
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theory.f] In the following we let the vectors n represent ordered integer pairs (n x , n y ) 
such that | nasi , \n y \ < A" max . The parameter A" max corresponds with a momentum 
cutoff scale of A = N max n/L. The modal field Hamiltonian has the formQ 



U = V r_I_fl 9_ ,l( ^ _ ti\ 6b_ g_ , _J8_ (s\ 2 1 / fi ooN 

~~ [ 29 1-adq n 2 \^ L 2 2 y (2L) 2 4! ^ (2L) 4 "™J l^U.ZO; 
n 

~l~ (2L) 2 ff ^ ] Q.n 1 < ln2 a n 3 , ( lfi i 



ni +™2 +"3 +™4 =0 



where 



and 



6 = Ei' + (6-29) 



«n = V 3 7^ — T Z V- (6.30) 

ni,n 2 



In Figure |6.1j we have plotted the lowest energy eigenstates in the rest frame 
for iV max = 10 and L = 3tt (in units where \i — 1). This choice of parameters 
corresponds with 441 different momentum modes and a momentum cutoff scale of 



A = 3.33yU. The states shown in Figure |TT] are the three lowest eigenstates in the 



even and odd — > —<p symmetry sectors, and the energies are measured relative to 
the ground state energy. In our calculation QSE diagonalization was used keeping 
500 Fock states, and the stochastic error correction was computed to first-order using 
the series method. Error bars shown include statistical error and an estimate of the 
contribution from higher order corrections. We see clear evidence of a second-order 
phase transition near % = 0.9.0 We have labelled the energies of the states according 
to their physical interpretation in the symmetric phase. E\ is the energy for the one- 
particle state, E 2 (E 3 ) is for the two(three)-particle threshold, and E' 2 (E' 3 ) is for the 
first state above the two(three)-particle threshold. At finite volume these energies 
are continuous functions of the coupling g. One feature which was also observed in 



3 We refer the reader to j£J for a short introduction. 
4 Counterterms were calculated using finite volume perturbation theory. 
5 A more complete discussion of the critical coupling as well as critical exponents 
can be found in |J2| . 
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Figure 6.1: Energy eigenvalues as functions of jj as calculated by QSE diagonalization 
with first-order error corrections. 



^i+i HH * s ^ ne crossing of energy levels due to the double degeneracy of states in the 
broken symmetry phase. E 3 is connected to a one-particle state in the broken phase 
while E' 2 is connected to a two-particle state. The levels E' 2 and E% therefore cross 
near the critical point. 

Another interesting phenomenon is the appearance of a bound state in the broken 
symmetry phase. In both the odd and even symmetry sectors we can measure the 
ratio of the two-particle to one-particle energies: 

Table 1 



9 
4! 


E' 2 /E 2 


E' 3 /E 3 


0.2 


2.01(4) 


1.98(4) 


0.3 


2.01(4) 


2.05(4) 


0.4 


1.95(4) 


1.96(4) 


0.5 


1.87(4) 


1.87(4) 


0.6 


1.86(4) 


1.82(4) 



These results are consistent with the binding energies reported in 
indicate a ratio of 1.83(3) near the critical point. 



and [[|, which 
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6.5 Compact U(l) in 2+1 dimensions 



Compact U(l) lattice gauge theory in 2 + 1 dimensions is a simple but phenomeno- 
logically interesting gauge model. It is asymptotically free and in the usual continuum 
limit describes massless non-interacting photons. On the other hand if the contin- 
uum limit is reached by rescaling the mass gap to remain constant, one instead finds 
a confining theory of massive bosons [111. The Hamiltonian has the form 



h = - ~ 2x J2 cosdp - ( 6 - 31 ) 

i 1 p 

In (|6.31|) A[ are link gauge fields, 6 P is the sum of the links circuiting a plaquette, 

9 P = A h + A l2 - A h - A u , (6.32) 



and 



x 



-4 -2 

e a 



(6.33) 



is the strong coupling parameter, which tends to infinity as the lattice spacing a 
goes to 0. We follow the notation of |H] in which an overall constant factor of y 
multiplying the right-hand side of ( |6.31| ) is suppressed. The energy levels we measure 



are therefore in units of y. 



The diagonalization of lattice gauge Hamiltonians is constrained by the require- 
ments of gauge invariance. To preserve gauge invariance it is most convenient to use 
a basis which diagonalizes the electric field part of the Hamiltonian 



Ed 2 

i 



\n v ) 



M ■ 



(6.34) 



As our next example of stochastic error correction we will address the 4x4 lattice 
system at x — 1 using this electric field basis. In it was noted that this poses 
a challenge to standard diagonalization techniques. Even on the small 4x4 lattice 
a surprisingly large number of states, about 10 7 ~ 10 8 , are needed to accurately 
describe the low energy spectrum at x — 1. This problem can be circumvented by 
modifying the basis states to incorporate more of the physics of the ground state. For 
example one can introduce a disordered background of magnetic flux as suggested in 
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22 1, and that approach is followed in an ongoing project |37| . However we would 
like to directly address the problem described in [215 and show how the stochastic 
Lanczos method handles the proliferation of large numbers of basis states in the 
original electric field basis. 

We will choose our starting subspace S to include all basis states 



\n v ) 



(6.35) 



which satisfy 



< 



(6.36) 



and which can be reached from the strong coupling vacuum by at most two transitions 
via the plaquette operators exp(±z#p). We take the auxiliary space Z to be the 
subspace spanned by basis vectors 



2 

max' 



(6.37) 



Using matrix diffusion Monte Carlo, we diagonalize the subspace formed by the states 



\d,j) = [P z (H-a)P z ) d \j). 



(6.38) 



\j) are the eigenvectors of the Hamiltonian restricted to the original subspace S. In 
our calculations we use a - 12 



L 2 max and d 



0, 1, • • • 12 for cutoff values, 



24,28,32. 



(6.39) 



In Table 2 we show the results for the ground state energy Eq for the different cutoff 
values i^ax an d the extrapolated value at L^ax = 00 • The errors shown are estimated 
statistical errors. For comparison we show the results of |21|] obtained using Green's 
function Monte Carlo (GFMC). 



Tabic 2 





L 2 — 24 

"max '"^ 


^max — 28 


^max — 32 


^max — 00 


GFMC 


E 


-7.394(3) 


-7.430(3) 


-7.438(3) 


-7.442(4) 


-7.4432(5) 



In Table 3 we show the masses for the lightest six particles in the system extrapolated 
oo. We have labelled the particles according to their spin J and 



to the limit -Z^^ax 
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sign under conjugation C : A — > —A. We also include results from j2l] for the lowest 
antisymmetric and symmetric glueballs. 

Table 3 



J c 


Mass 


GFMC 


|0"> 


3.03(2) 


3.01(6) 


|0+> 


4.03(3) 


4.05(8) 


|2"> 


6.8(1) 




|2+> 


6.8(1) 




|0+> 


7.0(2) 




|0"> 


7.1(2) 





The results we find appear in agreement with [21|. Unlike most Monte Carlo al- 
gorithms, the SEC method is able to find excited states with the same quantum 
numbers as lower lying states. This was also evident in the </>2+i example where we 
could track many different states crossing the phase transition. The reason for this 
advantage goes back to the design of stochastic error correction as a Monte Carlo 
improvement of a diagonalization scheme. For the U(l) example one can reliably 
find the eigenvalues and eigenvectors for the first twenty or so states in the low energy 
spectrum. 



6.6 Hubbard Model 

The last example we consider is the two-dimensional Hubbard model defined by 
the Hamiltonian 

H = - t E ( c - c ^ + 4o + u ZH^W)- ( 6 - 4 °) 

<i,j>; ar=1,i i 

The summation < i,j > is over nearest neighbor pairs. c\ a (ci a ) is the creation(annihilation) 
operator for a spin a electron at site i. t is the hopping parameter, and U controls 
the on-site Coulomb repulsion. The model has attracted considerable attention in 
recent years due to its possible connection to (i-wave pairing and stripe correlations 
in high-T c cuprate superconductors. In spite of its simple form, the computational 
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difficulties associated with finding the ground state of the model are substantial even 
for small systems. Fermion sign problems render Monte Carlo simulations ineffective 
for U positive and away from half-filling, and the collective effect of very large num- 
bers of basis Fock states make most diagonalization approaches very difficult. A brief 
overview of the history and literature pertaining to numerical aspects of the Hubbard 
model can be found in []53] . 

In terms of momentum space variables, the Hubbard Hamiltonian on an iV x iV 
periodic lattice has the form 

H = - 2t E ( cos ^ + cos ^) [ c lU c L Py + 4l Pv ,A*, Py ] ( 6 - 41 ) 

Px,Py=0,— ,N—1 

~ N 2 / j Px,P y gx,g y r x ,r y ^s x ,s y - 

Px—qx+r x —s x =0 mod N 
Py~1y J <- r y~ s y = Q mod N 

As a test of our methods, we use QSE diagonalization with stochastic error correction 
to find the ground state energy of the 4x4 Hubbard model with 5 electrons per 
spin. The corresponding Hilbert space has about 2 • 10 7 dimensions. For the QSE 
diagonalization we use momentum Fock states which diagonalize the quadratic part of 
the Hamiltonian. The Hamiltonian is invariant under the symmetry group generated 
by reflections about the x and y axes, interchanges between x and y, and interchanges 
between J, and |. We find it convenient to work with symmetrized Fock states. We 
will compute stochastic error corrections to first order using the series method. 

In Table 4 we present results for the ground state energy. We encountered no 
trouble with the sign problem, and in fact one can easily see that each term in the first 
order series expression (|6.5|) is negative definite. The energies are measured relative 
to the energy of the Fermi sea at U = 0. The errors reported are statistical errors 
associated with the first order SEC calculation. Where available, we compare with 
the results presented in [^5], which we label as Exact, Projector Quantum Monte- 



Carlo (PQMC), and Stochastic Diagonalization (SD). Stochastic diagonalization is 
a subspace diagonalization technique similar to QSE but one which uses a different 
method for selecting the subspace and is based on a variational principle Al- 



though the precise number of basis states used in the SD calculations is not listed, 



Chapter 6: Stochastic error correction 



81 



we infer from numbers reported for a modified 4x4 Hubbard system that roughly 
10 5 states were used.Q 



Table 4 



Coupling 


States 


QSE 


QSE+SEC 


Exact 


SD 


PQMC 




100 


-.4797 


-.50147(5) 








U = 2t 


500 


-.4945 


-.50181(3) 


-.50194 


-.5010 


-.44(5) 




1000 


-.5006 


-.50198(1) 










100 


-1.620 


-1.8113(4) 








U = 4t 


500 


-1.748 


-1.8242(3) 


-1.8309 


-1.829 


-1.8(2) 




1000 


-1.800 


-1.8302(1) 










500 


-2.558 


-2.7073(4) 








U = 5t 


1000 


-2.651 


-2.7208(2) 


-2.7245 


-2.723 


-2.9(3) 




2000 


-2.685 


-2.7231(1) 









Apparently QSE diagonalization with SEC handles the 4x4 system quite well with 
relatively few states. Much larger systems are being studied using both higher series 
corrections and stochastic Lanczos techniques 150 . 



6.7 Summary and comments 

In this paper we presented two versions of stochastic error correction, the series 
method and the stochastic Lanczos method. The series method starts with eigen- 
vectors of the Hamiltonian restricted to some optimized subspace and includes the 
contribution of the remaining basis states as an ordered expansion. The stochastic 
Lanczos method starts with eigenvectors of a Hamiltonian submatrix and constructs 
matrix elements of Krylov vectors using matrix diffusion Monte Carlo. This method 
has the advantage that the starting vectors need not be close to exact eigenvectors. 

We presented three different examples which demonstrate the potential of the 
new approach for strongly coupled scalar, gauge, and fermionic theories. In the first 

6 The discrete symmetries of the system were not utilized in their calculations. 
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example we calculated the low energy spectrum of 4>2+i using the series method, and 
in the second example we found the spectrum of compact U(l) in 2 + 1 dimensions 
using the stochastic Lanczos method. In both examples we found agreement with 
results from the literature. We also found that unlike typical Monte Carlo results, 
the SEC method is able to find the eigenvalues and eigenvectors for excited states 
with the same quantum numbers as lower lying states. This advantage is due to 
its design as a Monte Carlo improvement of a diagonalization scheme. In the last 
example we found the ground state of the 2+1 dimensional Hubbard model using QSE 
diagonalization and first order series stochastic error correction. In this calculation 
we encountered no fermion sign problem and found that our methods yielded very 
accurate results with far less effort than existing techniques. We believe that the 
methods we have presented hold considerable potential for studying a wide range of 
non-perturbative quantum systems and answering questions difficult to address using 
other methods. 



Chapter 7 



EQSE Diagonalization of the 
Hubbard Modeli 

7.1 Introduction 

The enhanced quasi-sparse eigenvector (EQSE) method of solving quantum field 
theory Hamiltonians is the combination of quasi-sparse eigenvector (QSE) method 
|39| with a stochastic calculation for the contribution of the remaining basis states 
§40|| . The Hubbard model was chosen as a laboratory for testing this method for 
several reasons. First, we believed the approach could yield results. The basis vectors 
can be specified in a few words of data and the Hamiltonian is sparse is momentum 
space. On the other hand, the ground state of the Hubbard model is known for its 
inclusion of an extraordinary number of Fock states JU so the model presents a non- 
trivial challenge to the quasi-sparse approach. Finally, the Hubbard model is thought 
to be a physically relevant model for superconductivity [I]]. While the description of 
the model is simple, solutions have been difficult and any promising new approach is 
worthwhile. 

^ucl. Phys. B (Proc. Suppl.) 90 (2000) 202-204 
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We work on a 2- dimensional spatial lattice with the Hubbard Hamiltonian 

H=-t ( c l c i- + sW) + u J2( c l^4i c n) 

<X=U 

< i, j > nearest 

There are 2 species of electron so there are 4 possible states for each lattice site. 
Thus the 8x8 Hubbard model has 4 64 dimensions. Even after using particle conser- 
vation to partition the space there are more than 10 32 basis vectors. In this large 
space the Hamiltonian is clearly sparse But the equality of the off-diagonal elements 
contributes to an extraordinary number of Fock states in the ground. 

It is thought that the D-wave correlator Cd 2 _ 2 (r) is an important indicator of 
superconductivity |25[ . 



7.2 Hamiltonian Momentum Lattice Formulation 

After making space periodic we use the Fock states as our basis set. The Hamil- 
tonian conserves momentum which helps limit the number of relevant basis states. 

7171 7T P\ t 

#kin = ~2t 2^( COS ~T + COS -^) a np,a a np,a 

k — I + m — rt = 
p — q + r — s = 

There is a 16 member symmetry group generated by reflections in the x and y planes, 
x <-> y andJ,-*-»|. For the purpose of finding the ground state, we use only symmetrized 
basis states. 

The first step in the calculation consists of picking a basis set of size iV and 
diagonalizing its submatrix of the Hamiltonian using the Lanczos method [ffE. The 
ybasis vectors which least contribute to the ground are then discarded, replacements 
are chosen and the Hamiltonian is again diagonalized. When the ground energy Eq 
obtained in this manner converges the QSE step is complete. 
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The EQSE is a first order correction to this result. We calculate it stochastically. 
Let C be a set of basis vectors for the complement of the iV-dimensional subspace. 
Then, choosing representative basis vectors v G C with probability P(v) . 



E = E + Average \ ; , - 1 



where |0) is the ground state of the N dimensional subspace. The expectation is thus 

FO(v) 

,-_( • 

7.3 Results 



Results were obtained for the ground state energy, wavefunction and <i-wave cor- 
relator. The computing time was about 2 days on a 350Mhz Pentium II. Where avail- 
able, we compare with Husslein et al [25| results labeled Exact, Projector Quantum 



Monte-Carlo (PQMC), and Stochastic Diagonalization (SD) (which uses a different 
method of choosing the subspace than QSE). 
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Ground State Energy 4x4 Hubbard Model occupied) 



( jOHTll 111 (X 


kJ uduCQ 


Vc/kJ-LJ 


F,OSF, 

J_JV/ kJ-LJ 


UAdL L> 


SD 






oU 


-.4/4/ 1 


-.oUlz / (0 J 








T T O 

U=2 


100 


-.47967 


-.50147(5) 


-.50194 


rni 

-.501 


a n 

-.49 




500 


-.49454 


-.501ol(o J 










1UUU 


-.oUUoz 


-.oUlyo 










50 


-1.5f0f 


-I.o0oo5(5 ) 








TT A 

U=4 


1UU 


-l.ozUo 


-l.ollo(4 ) 


-l.ooUy 


1 son 


-1.8(2) 




500 


-1.7476 


-1.8242(3) 










1000 


-i.800o 


-1.8302 










50 


-2.2450 


-2.6663(8) 








T T r 

U=5 


100 


-2.6622 


-2.6724(7) 


-2.7245 


O TOO 

-2.723 


-2.9(3) 




500 


-2.55 10 


-2. / ( o(4 J 










1000 


-2.6512 


-2.7208(2) 










2000 


-2.685 


-2.7231 










50 


-2.963 


-3.615 








U=6 


100 


-3.103 


-3.635 










1000 


-3.452 


-3.697 










2000 


-3.595 


-3.723 
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Ground Energy 8x8 Hubbard Model(g| occupied) 



( nnn lTin 1 

v^uupiiiig, 


Stcitcs 








5U 


-.281 


-2.58(2) 




1UU 


A AO 

-.443 


O A m Zo\ 

-2.49(2) 


T T O 

U=2 


500 


-1.221 


O /j M / O \ 

-2.40(2) 




1UUU 


-1. ( 01 


-z.4Ub(o ) 




2000 


1 Mr/ 1 

-1.95b 


O 1 O O / M \ 

-2.423(9) 




4000 


1 m r o 

-1.958 


O A OT 

-2.427 




r - m 

50 


O 1 1 

-.811 


m 1 o / n \ 

-9.18(b) 




i m m 

100 


1 o o n 

-1.38b 


-8.48(b) 


T T A 

U=4 


r mm 

500 


O A A M 

-0.449 


-7.58(5) 




i mmm 

1000 


A ^7M O 

-4.798 


-7.59(3) 




OM MM 

2000 


-5.374 


-7.620(4) 




/) MMM 

4000 


-5.387 


-7.621(5) 




r m 

50 


1 OO o 

-1.222 


-13.33(8) 




1 MM 

100 


-1.88 


1 o oo/o\ 

-12.33(8) 


T T r 

U=5 


r - mm 

500 


-4.b5 


-10.65(6) 




1 MMM 

1000 


f A A 

-b.44 


a m ^r/o\ 

-10.65(3) 




OMMM 

2000 


t oo r 

-7.225 


-10.671(5) 




4UUU 


- ( .z3b 


-lU.bbz(o ) 




50 


-1.6 


-18.2 




100 


-2.4 


-lb.b 


U=6 


500 


-5.9 


-13.93(8) 




1000 


-8.15 


-13.78(4) 




2000 


-9.110 


-13.888(7) 




4000 


-9.113 


-13.880(8) 



The d x 2_ y 2 correlator was also obtained using the QSE algorithm. The 4x4 re- 



sult again matched that of Husslein et al pqj . The EQSE calculation has not been 
completed and we therefore omit the data. 
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7.4 Conclusion 

As we can see, the ground state for the 4x4 Hubbard model can be well described 
with about 1000 symmetrized states and the 50 state results yield remarkable accuracy 
when the first order correction is included. In the 8x8 case it is clear that even 4000 
states are not sufficient to describe the ground state. The precision of the first-order 
values will be determined with the completion of higher order calculations. 

Further advances will come in the refinement of the enhancement technique. Bet- 
ter importance sampling will speed convergence of second and higher orders contri- 
butions. Other extensions will be calculation of excited states of the Hamiltonian, 
correlation functions, binding energies and other quantities of interest. 



Chapter 8 

Higher order calculations in the 
Hubbard model 

8.1 Introduction 

The last chapter contained an application of the quasi-sparse eigenvector approach 
to determining the energy of the Hubbard model. The work was done with a Fock 
state basis. The 4x4 Hubbard model has few enough states that it is possible to 
directly diagonalize the Hamiltonian and we were able to compare our results with 
the exact answer. 

The accuracy of the results obviously depended on the dimensionality of the ex- 
actly solved subspace, Q. We summarize here the results for a U/t ratio of 4, which 
may be in the physically relevant range for High-Tc superconductivity. A subspace 
spanned by 50 symmetrized states yielded the ground state energy to within 15%. 
This slowly declined as we increased the dimensionality. A dimensionality of 1000 
was sufficient to within 2%. On the other hand, a first order perturbative calculation 
(EQSE), when applied to the 50 state space was also sufficient to achieve 1% accuracy 
and when applied to the 1000 state space the answer was accurate to 4 significant 
figures. 

For the 8x8 Hubbard model the situation was much less clear. While quasi-sparse 
and first-order corrected results were obtained, their accuracy could not assessed. 



89 



Chapter 8: Higher order Hubbard 



90 



The EQSE results appear to converge for large initial subspaces. However, a close 
look shows that the uncorrected quasi-sparse result also appears to converge for Q 
over about 1200 states. The change from 2000 to 4000 states only lowers the ground 
energy by about .01 units. But the first order result is more than 2 units lower still. 
This is evidence of the extreme number of Fock states in the Hubbard ground. Even 
if there is no further falloff in contribution, one could estimate from this that a quasi- 
sparse calculation with 400,000 symmetrized basis states in Q would be required to 
reach the accuracy of the the first order result. 

The purpose of this paper is to further explore the ground state energy of the 
Hubbard model using several approaches including variations of the ordered expansion 
and the stochastic Lanczos method presented in chapter 6. 



8.2 Calculations 



On a 4x4 lattice with 5 occupied sites for both spins, the perturbative ground is 
given by 



2 . . . . 
1 

®<g)<g). 
-1 . <g . . 
-10 12 
t modes 
A typical excited state is 

2 . . . . 2 

1(g)®.. 1 

<8 <8 . . 

-1 . <g) . . -1 

-10 12 

t modes 



2 . . . . 

1 .(g).. 
<g> <g> <g> . 
-1 . <g> . . 
-10 12 
| modes 



|0> 



a-iialoa-i-ialo|0) 



-10 12 
| modes 

We can define a distance function of a state as the number of differences from the 
perturbative ground. The example above has 2 differences among the j modes and 2 
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differences among the I modes. We say its distance is (2,2). In cases where we are 
only interested in the greater of the up and down distances we would say the distance 
is 2. 

When the ground state energy computed in chapter 7 was considered as a function 
of the number of states, we found that the slope would suddenly drop off at around 
1200 states for the 8x8 Hubbard model. For the 4x4 case, the drop in slope occurred 
at around 30 states. The fillings we have used, 5/16 in the 4x4 case and 25/64 in 
the 8x8 case both have a non-degenerate U — ground state. On close examination, 
it became clear that the first states found in the quasi-sparse diagonalization were 
always distance (2, 2). When these states, for which H has non-zero matrix elements 
with the free ground state, were exhausted, the slope of the energy curve would flatten 
out. 

The calculations presented in this chapter all used the natural cutoff of distance 
< 2 to define Q, the base subspace. This simple definition saved time in determining 
if new states were contained in the base subspace. Further, given this definition of Q, 
a clear relationship between powers of H and distance of a state became clear. H has 
non-zero matrix elements only between states whose relative distance is less than or 
equal to 2. As an example, in a calculation of H 4 between the base space and itself 
the first and third intermediate states would be states with maximum distance 4 and 
the second intermediate state would have maximum distance 6. 

The first extension of the previous results was inclusion of higher order contribu- 
tions in the ordered expansion. This was done in two ways. Second order Monte Carlo 
calculations were attempted using the ordered expansion described in chapter 6. The 
results were on the same order as the first order results and are not presented here. 
Higher order calculations in this expansion were difficult because of the many dif- 
ferent terms that contribute. Instead, the Brillouin-Wigner degenerate perturbation 
method was used to calculate higher terms. This approach yielded some interesting 
results which will be discussed. 

The other method of extracting data was to use Monte Carlo to calculate the 
matrix elements of powers of the Hamiltonian. Once these are known there are 
several ways they can be used. The stochastic Lanczos method uses these matrix 
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elements to generate a matrix for if in a Krylov subspace generated from the ground 
of the base subspace. In theory one could diagonalize the Krylov subspace to obtain 
the eigenvectors and eigenvalues of H. The usefulness of this method was limited by 
the very small determinants of the matrix thus created and the uncertainties of the 
matrix elements calculated using Monte Carlo. 

A more careful choice of subspace, using states generated by higher eigenvectors 
of Q may solve this problem. Although some work was done in this direction, it is not 
presented here. Instead, we used the power method to extract the eigenvalue using 
the same matrix elements for H. 



8.3 Brillouin-Wigner perturbation theory 

The Brillouin-Wigner perturbation method is described in Appendix A. The gen- 
eral idea is to solve a self-consistent equation for the eigenvalue. For the desired 
eigenvector \ip) in Q, we find the eigenvalue and components as follows. 



Hij = (i\H\j) 

A (i) _ v (t|H|fc)(fc|H|t) 

^ij ' l^\k)iQ E-{k\H\k) 

A2) v (i\H\k)(k\H\l)(l\H\i) 

V ^ \k)Al)tQ (E-(k\H\k))(E-(l\H\l)) 

\k) + 10 

A (3) _ v (i\H\k)(k\H\l)(l\H\m)(m\H\i) 
"ij - \ k ),\l),\m)iQ {E-{k\H\k))(E-{l\H\l))(E-{m\H\m)) 

\k)^\l),\l)^\m) 

etc. 

This could be solved by feeding the n-th order result for E into the (n + l)-th 
order calculation's denominators. Alternatively we could keep a running estimate for 
E as the Monte Carlo proceeds and, at each step, use the current estimate. Instead, 
we chose to calculate the left hand side E t as a function of the right hand side E r . 
We then take the intersection of the graph of this function and the line Ei = E r as 
the eigenvalue. 
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Each Monte Carlo trial consists of a trajectory beginning and ending in Q. At 
each step the weight is adjusted by the appropriate probabilities and, outside of Q the 



appropriate denominator is divided. Figure |8.1| presents data for the 6x6 Hubbard 
model with E r = 7. We would have to run this again with E r = 5 to draw a curve 
and find the actual eigenvalue. Two of the runs used a distance cutoff of 6 and two 
used a distance cutoff of 8. 

A guiding scheme was used to increase the likelihood that paths return to Q and 
to prevent them from going beyond the cutoff distance. Were it not for the denomi- 
nators the guiding would be perfect in the sense that all paths would have an equal 
contribution to the matrix elements. A better guiding scheme would take account of 
the diagonal energy and better sample those states with smaller denominators. 

The convergence is excellent up to order 4 and is good up to order 6 or 7. After 
that we consistently find values for the minimum eigenvalue that are much too low. 
One can account for this by the non-linear process of diagonalizing that is used to 
find the eigenvalue from the Aij matrix. The distance cutoff 6 and distance cutoff 8 
results remain close at order 5 and 6. After order 7 the noise in the cutoff=8 results 
becomes dominant. 



8.4 Power method 

The power method is very simple in theory. We start with the ground state of H 
restricted to Q, run trajectories with H — a and tabulate the results. If the energy 
shift, a, must be chosen such that \E — a\ > \Ej — a\ for all other eigenvalues Ej, 
then the ratio of successive powers of H will approach the value (Eq — a). 

It is not surprising that the calculation suffers from a sign problem. The guiding, 
described in the previous section, results in the contributions of different paths at 
each order being close to each other in magnitude but of varying sign. Functionally, 
after around 7 or 8 steps this problem significantly adds to the convergence time. 

Several modifications were made to minimize the impact of long paths. The first 
was that only paths outside of Q were tabulated. By restricting to such "connected" 
paths we decrease the contribution of long paths. Let Q be the projector onto the 




Figure 8.1: Brillouin-Wigner calculation for energy as a function of order. gnu6_2 
and gnu6_3 have max distance 6. The others have max distance 8. 
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space Q and P = 1 — Q. And let C n ^ be the n-step matrix element between i and j 
that we get by restricting to paths outside of Q. 

Ci = QHQ 

C n = {QHp){pHP) n - 2 {PHQ) 

H n = C n + Cn^H 1 + ... C x H n - x 

The disadvantage of this approach is that it increases the memory needs. 

If an intermediate state Q is chosen with known probability we can calculate 
the contribution of paths that pass through it at the m-th step. We multiply the 
Monte Carlo results for (i\H m \ip) and (ip\H n \j) together to find the contribution to 
H™ +n . Effectively, this allows the calculation of twice as many steps before the sign 
problem kicks in. The concern with this approach is that more time is spent on each 
state at the m-th step so the space of possible m-th states is not as well sampled. If 
the contribution to C m+n is a relatively smooth function of the m-th state this is not 
a problem. In practice, this method was only worthwhile for higher order terms. 



We present the results for the 8x8 Hubbard model in figure |8.2| . The distance 
function was limited to a maximum of 8. We plot the ratio of H n /H n ~ 1 for n — 1 to 
11. The energy shift is set to 30. Two trials were run. C\ through Cj were calculated 
with paths starting in Q and Cg through Cn used the intermediate state method. 
As expected, the paths approach a decaying exponential. Using gnuplot to fit points 
5 through 11 yields an asymptotic value of —39.74 for one trial and —39.75 for the 
other. Other choices of points to fit result in small changes in these results. Overall, 
for a distance cutoff of 8, the 8x8 Hubbard model yields a ground state energy of 
-9.74 ± .05. 



We present results for the 6x6 Hubbard model in figure |3.3| . With a distance 
cutoff of 6, both curves are on top of each other and yield an energy of -4.887. With 
a distance cutoff of 8 the answer is -4.94. In these calculations the intermediate state 
method was not used so the results at n > 7 show larger errors and are not included 
in the fit. 
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8.5 Conclusions 

For the computing power we had available, the 8x8 Hubbard model poses a serious 
challenge. We have been able to extract a good estimate for the ground state energy 
in the || filled sector. The same techniques may also work for calculating correlators 
in the ground state. The most powerful approach we have used is the power method 
but the matrices C n extracted for this method are also available for other methods of 
analysis such as the modified stochastic Lanczos method alluded to above. While the 
Brillouin-Wigner method may hold promise for the future, it was less precise for our 
purposes than the power method. The non-linearity of the inversion step makes it 
particularly sensitive to noisy data. Improvements in guiding or sampling strategies 
and filters to allow linear treatment of noisy data may help both the Brillouin-Wigner 
and power method approaches. 



8x8 Hubbard model, U=4 



method 


QSE first order correction power method 


states 


2000 


4000 


2000 


4000 


1200 


E 


-5.374 


-5.387 


-7.620 


-7.621(5) 


-9.74(5) 



A brief look at the summary of results for the ground energy shows that despite 
the apparent convergence of the first-order correction, the answer is still far from 
correct. Computations for distance cutoffs of 6 and 10 will allow a determination of 
the precision of the power method result presented here. 

Another interesting extension will be the calculation of ground state energy for 
other fillings. Exploring the fillings near the one presented here will allow determi- 
nation of the energy to add or remove an electron from the lattice. 
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1 r 

fa(x) 

fb(x) — - 
fc(x) 
fd(x) 
'R6x6d6A' + 
'R6x6d6B' 
'R6x6d8A' * 
'R6x6d8B' :---B 




-34 - 



-34.2 - 
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10 



12 



14 



Figure 8.3: 6x6 Hubbard model: Ratio of (ip\H n \ip) over {ip\H n ~ 1 \ip) as a function of 
n. Two trials with max distance 6 and two trials with max distance 8 are shown as 
well as their best fit curves. 
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Appendix A 

Degenerate Perturbation Theory 



Let |1), |2) ... \n) . . .be eigenstates of H Q with energies E n . Assume E 1 through 
E m are near each other but are all sufficiently far from E m+1 , E m+2 etc. And let 
be an eigenstate of (H + V) with eigenvalue W where W is close to the E\ through 
E m . 

By the eigenvalue equation 

(H + V)\if>) = W\i/>) 
which if we expand in terms of \n) becomes 

J2(H + V)\n){n\i/>) = Wj^ WW) 

n n 

So, taking a contraction with (i\ on the left we get 

Ei(i\il>) + ^2(i\V\n){n\7p) = W(i\^) (A.l) 

n 

This is just an (infinite dimensional) eigenvalue matrix equation. 

' E 1 + (1\V\1) (1\V\2) (1\V\3) 
(2\V\1) E 2 + {2\V\2) (2\V\3) 
(3\V\1) (3\V\2) E 3 + (2\V\3) 

(A.2) 
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= W 


m 
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If we temporarily fix the value of W this is just a homogeneous first-order equation 
in (n\ip). Eqn ( |A.2| ) will only have a solution for special values of W, but we can use 
it to determine (m + (m + 2\ip), etc. in terms of (1|^>) through {m\ip). 

If we treat as fixed for i < m we get the equation 

Ei<m(rn + l\V\i){i\iP) 
Zi< m (™ + 2\V\i)m 







(m + 1 -0) 






[S] 


(m + 2\ijj) 




where 










E m+1 + (m + l\V\m + 1) - 1 


[S] = 


(m + 2\V\m + 1) 





(m + l\V\m + 2) 
E m+2 + (m + 2\V\m + 2) - W 



Since W is not near E n for n > m + 1 and V is small, [S] must have non-zero 
determinant and be invertible. Since the right hand side is linear in for i < m 
we can ( by multipling by [S^ 1 ]) express (n\ip) for n > m + 1 as a linear combination 
of (#> • 

For n > m + 1 we therefore write 



(n\iP) = J2 A mi^ 

Plugging this into (|A.1| ) we get 
(W-E n -{n\V\n))J2 A m(W) = E^l^l^^l^)" 



i<m i<m 



i<m 



Y,^ n \ v \ n ') A ^ 1 

n' 7^ n 
n' > m + 1 



which can be solved for A ni by fixing coefficients of (i 

1 



Ani W-E n - (n\V\n) 



\{n\V\i)+ {n\V\n')A n ,)j 

n' 7^ n 
n' > m + 1 



(A.3) 
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Although this equation has A n n on the right hand side it could be solved iteratively 
since the A n n comes with an extra factor of V relative to the A ni on the left hand 
side. However, remember that A n i is really a function of W so we are not ready to 
solve it until we find W. 

Now we are ready to turn to the eigenvalue problem in the |1) through \m) sub- 
space. Once again using eqn ( |A.1| ) 

W(i\i;) = (£ i+ (i|V|i))<#)+( J2 (AV\j)(M))+ £ J>|V|nK;<jiV>> 

j ^ i n > m + 1 

j <m 



Ei + (t\V\t) + Z (z\V\n)A ni )(i\i;) 

n > m + 1 / 

+ E ... (<«W>+£ ^ (i\V\n)A nj )(M ( A - 4 ) 

j f i V n > m + 1 J 

j < m 

Equation (|A.4|) is an eigenvalue equation, ([M] — W r )[(z|'0)] = 0, for W with the 
m x m matrix [M\ given by 

M u = {E, + (i\V\i) + Yl (i\V\n)A n ^J 

n > m + 1 

and 



My 



n > m + 1 

The only problem is that [M] is itself a function of W through A ni . Again, however, 
this can be solved iteratively since the A n i terms all come in with an extra factor of 
V. 

So, the general procedure is to solve the eigenvalue equation to first order in V 
(ignoring all the terms with A ni since A ni itself has order 1) and then to use this W 
to solve for A ni to first order. (We can ignore the terms with A n >i on the right hand 
side since they are second order in V^.) Now we can continue by substituting the first 
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order value of A ni into [M] and now calculating W to second order. If we substitute 
this and the first order values for A n n into ( A.3 ) we can solve for A ni to second order. 

At order k we solve the eigenvalue equation for W using A n i to order (k — 1) and 
then use W to order k and A^ to order (k — 1) to solve for A ni to order k. Of course, 
there are m solutions for W at order k and we have to match the solution used for 
order (k — 1). 



